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 () 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, models are also scalable. This paper shows how 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 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.
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 () [4] and Continuous Petri Nets () [5].
In addition, Hybrid Petri Nets (), that combine all features of and were proposed in [1, 5]. provide a flexible tool to integrate smoothly continuous and discrete states. Furthermore, coloured Petri nets () [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 and , coloured Hybrid Petri Nets have been proposed in [8–10] combining all the merits of and into one net class. The different notions of firing rates of facilitate the modelling of multitemporal scales, while the use of colours allows a parametrised modelling of large systems with repeated components. In addition, can be converted algorithmically into equivalent low-level via an unfolding process [11]. Indeed, 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 () dynamics [12–15] is an important biological pathway that has been recently widely studied. Intracellular 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 dynamics. So far, many models were introduced to analyse 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 model cannot be achieved by using a single modelling paradigm (i.e., deterministic or stochastic).
On one hand, the 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 utilises [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 behaviour. The resulting model is flexible since it permits adjusting the model parameters to allow experimenting with different simulation configurations. In our modelling approach, 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 dynamics as well as the utilisation of for modelling biological processes. Section 4 shows how can be used to construct and execute biological models illustrated by the case study. Simulation results of the 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 () 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, 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 standard were proposed, such as Stochastic Petri Nets () [4], Continuous Petri Nets () [5], and Hybrid Petri Nets () [1, 25].
extend by associating each transition with a random firing rate. The semantics of an 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 () were proposed by considering different transitions types and different arcs types [27, 28]. Moreover, in order to model continuous concentration changes of biological processes, were proposed in [1]. The semantics of a 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, were proposed in [1] and further augmented to Generalised () in [25].
The motivation of , compared to , is to support more timing features, including stochastic, discrete deterministic, and continuous deterministic firing of transitions [25]. 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 , 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.
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 and interchangeably unless it is explicitly stated.
2.2. Coloured Hybrid Petri Nets
Coloured Petri Nets () extend low-level with a set of finite colour sets [6, 7]. offer a parameterised approach to construct large biological systems; there are many coloured net classes, among them coloured continuous Petri nets () [7] and coloured stochastic Petri nets () [31].
By augmenting with colours, we obtain a new class of hybrid Petri nets called Coloured Generalised Hybrid Petri Nets () [8]. In , 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.
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 to model diffusion in 3D space.
(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 Figures 1(a) and 1(b). Figure 1(a) illustrates a component at index and its six neighbours at different indices. These six neighbours are calculated by or , while Figure 1(b) describes a component at the index and its 26 neighbours at different indices. These 26 neighbours are calculated by or . 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). 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 , but also employs their colour expressions. To present a component we use the graph structure of 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 , , and , 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, , representing a cube form. Moreover, we define six variables, of type , of type , and of type , as well as three-dimensional coordinates, to address every component in space. Each of these components has a set of neighbours which are defined by the function : 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 modelling the diffusion to all positions fulfilling the guard .
In our developed 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 of 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 permit to combine different timescales.
2.4. Intracellular Calcium Dynamics
This section provides an example of how 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 dynamics, followed by a detailed discussion of the model construction.
2.5. Calcium Release
Intracellular calcium is organised by the release of ions from the cell membrane or internal stores such as Endoplasmic Reticulum (ER). Releasing of from ER occurs via the opening of inositol-4,5-triphosphate receptor () 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 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 release from the ER to cytosol.
As it is detailed in [17], () channels consist of four identical and independent subunits which can be combined with the deterministic part to study the 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 , one for activating , and a third for inhibiting. Each of these sites is either busy or idle by their particular regulatory agonist ( or ) [17, 34]. Building on the experimental studies in [35], Marchant et al. observed that in the absence of bound , does not bind to its activating site. Therefore, the two states of unbound 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 (the fastest), binding by activation, and binding by inhibition (the slowest). As dissociation and association is faster than the competing processes, transitions between states of unbound and bound are the first to reach the equilibrium state. Therefore, the unbound state is grouped with the bound 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 state.
Finally, the opening of an channel is influenced by the neighbouring 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 dynamics in two-dimensional space, while in [20–23] diffusion is studied in three-dimensional space. Furthermore, a hybrid model was introduced in [16] to study the 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 an 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 an 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 . Figure 3 provides the complete model.
The construction of the 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 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 the 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 as shown in Figure 3, where the place gives the calcium concentration in ER, gives the calcium concentration in the cytosol, and the places and give the calcium concentration in mobile and stationary buffers, respectively.
Seven (coloured) continuous transitions represent the interactions between in cytosol and ER as well as buffers. Calcium inflow from the into occurs through the transitions and . On the contrary, calcium can return back to the ER from the cytosol by a transition called . The interaction between the buffers and cytosol is controlled by the transitions , , , and . Since our main objective is modelling and simulating cytosolic inside a cube, we associate all places with the colour set , which is a compound colour set of the dimension (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 , , , and (see Figure 3) in the three-dimensional space.
We use a similar approach to the one presented in Figure 1 to model the diffusion dynamics. Therefore, the resulting model is flexible as it permits to adjust the number of grid positions by changing the constants that define the colour set , namely, , , and . In other words, the spatial dimension of the model is defined and later adjusted via the constants , , and .
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 our model, channels can be grouped into clusters.
The dynamics of the stochastic part of our is based on the one given in [36]. In [36], the authors’ goal was to construct a simple mathematical model for that can effectively describe the properties of release events by assuming that every has four equal and independent subunits. Each subunit has three sites of binding: a activation binding site, an 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, denotes the state of a subunit, where i characterises the state of the binding site ( means not bound), j the state of the activation binding site, and k the state of inhibiting binding site. However, according to [36], the states and 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, where denotes a state in which all binding sites are on, while the place 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 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 the colour set and the tuple , providing the channel index in 3D. Similarly, for modelling a cluster in a 3D grid, we use the colour set and the tuple as cluster index. However, channels are scattered in clusters. Therefore, we use the compound colour set 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 .
3.3. Modelling the Effect of Channel States on Calcium Releasing
Releasing of 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 was created with colour set to decide the channel state (open/close). The read arc between the place and the immediate transition 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 . The inhibitor arc between the place and the immediate transition prevents the repeated firing of the transition, as soon as there is a token on the place . On the contrary, another immediate transition will switch off the channel state, when the place 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 . Firing of the transition increases the marking of the continuous place , but in a specific region. The determined grid region, which permits the flow of from ER into cytosol through the transition , 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 which is given by [15, 19, 20].where 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 is a Boolean expression having the form:where , , and are the spacing between two grid positions in the x, y, and z directions, , , and give the position of the first cluster in the x, y, and z directions, , , and give the spacing between two clusters in the x, y, and z directions, and 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 , , , , , , , , and 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 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 place was added with the colour set , for providing the total number of open channels inside every cluster. The transition can be considered as the controller for the channel opening where tokens are added when the transition fires. In the contrary, the transition switches off the channel as soon as it fires. Besides, the tokens in the place are subsequently used by the deterministic regime to regulate the flow.
3.4. Modelling the Effect of Calcium Releasing on Channel Transitions
The four transitions of channel subunits are affected by concentration in the cytosol close to the channel, as discussed in [16, 36]. Therefore, it is required to sum up all 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 to sum up for each channel the nearby 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 constants L, 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 constructed model shown in Figure 3 by performing different experiments; it shows how easily we can scale 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 examine 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 .
The 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. 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 model result is valid for the simulation of a single channel in a single cluster.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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 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 concentration. Please note that although there is no release due to the channel being closed, there is transport due to the leak term.
4.2. Calcium Diffusion
The dynamics of intracellular calcium is described by diffusion. Therefore, we have performed two experiments to study the diffusion dynamics of the model. The first experiment was configured with only one cluster including 24 channels by setting , , and , while the second experiment was configured with four clusters with every cluster containing 12 channels by setting , , and .
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 Figure 6. 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 concentration inside the grid (cube form) when six channels are in the open state, while Figure 6(c) shows the same results of the concentration when seven channels are in the open state. Figure 6(d) shows 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 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 starts increasing again (Figure 6(f)).
(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 Figure 7. 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).
(a)
(b)
(c)
(d)
4.3. Validation of the 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 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 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 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, 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 the 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 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 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 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 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 concentration inside the cube. A comparison of the behaviour of our model with results known from similar models in the literature is provided in Table 2.
4.4. Simulation Procedure
The 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 can be found in [42].
Simulation of in Snoopy is performed by first unfolding the model into the equivalent one, followed by choosing a suitable hybrid simulation algorithm to execute the unfolded , 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 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], 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 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 .
5. Conclusions
This paper has presented the construction and simulation of 3D models using . 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 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 is showed by changing only a few parameters to run more than one experiment. The results of the developed have been validated against results known from the literature.
Our 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 to other modelling scenarios in future work.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Supplementary Materials
A Graphical Approach for the Hybrid Simulation of Diffusion 3D Bio-models. (Supplementary Materials)