Deterministic and Stochastic Hybrid Simulation

stochastic model contains random mechanisms and, hence, the simulation results can change in any run performed with the same input parameters.

Hereafter, the considered stochastic models are stochastic Petri nets (SPNs) already introduced in Section 4.3. The implementation of the stochastic transition by means of the Modelica language is part of Section 6.1.3. The concept is also applicable to an xHPNbio if it contains at least one stochastic transition.

Each transition of a SPN represents a reaction of a biological process and the places are the biological compounds which interact with each other. The tokens quantify the amount of each compound, for example, the number of molecules and the arc weights represent the number of units which are consumed and produced, respectively, in a reaction, for example, the stoichiometric coefficients of a biochemical reaction (see Definition 5.1).

The time when a reaction occurs is a random quantity defined by the hazard function (see Definition 4.34)

, Eq. 5-15

whereby is the putative firing time of transition which models the reaction, and

is an exponentially distributed random variable specified by the hazard function . This hazard function can be defined more precisely if a SPN represents a biochemical network.

Two possible hazard functions can be applied which depend on reading the tokens as molecules or as levels of concentrations (Heiner et al. 2008).

The first, called stochastic mass action hazard function, is directly derived from the mass action kinetics of continuous reactions (Wilkinson 2006). The mass action kinetics has to be modified such that the substances can be represented discretely by molecules instead of continuous concentrations. This general modification process is shown exemplarily by the following types of biochemical reactions.

A biochemical reaction of first-order has the form : → ,

whereby the constant represents the hazard that a molecule of will undergo the reaction and there are molecules of substance . This leads to the hazard function

⋅ . Eq. 5-16

A second-order biochemical reaction can have the form

: → ,

whereby the constant represents the hazard that a molecule of and a molecule of react with each other and there are molecules of and molecules of . Hence, ⋅ different pairs are possible so that the hazard function is given by

⋅ ⋅ . Eq. 5-17

It is also possible that a second-order biochemical reaction has the following form : 2 → ,

whereby the constant represents the hazard that two molecules of substance react. But only 1 /2 pairs are possible so that the hazard function is given by

⋅ ⋅ 1

2 . Eq. 5-18

This theory can also be applied to higher-order reactions and results in the following general form of the stochastic mass action hazard function regarding that the reaction is modeled by a SPN

⋅ →

∈

, Eq. 5-19

whereby is the transition-specific stochastic rate constant, is the token number of the input place of transition , and → is the weight of the arc → .

The second one, called stochastic level hazard functions, regards the tokens as levels of concentrations as introduced in (Calder et al. 2006). The concentration of each substance is then transformed to an abstract level. Therefore, it is assumed that the maximum molar concentration is and the number of the highest level is so that the amount of levels is

1. Then the abstract levels 0, 1, … , represent the concentration intervals 0, 0, 1 ⋅ , 1 ⋅ , 2 ⋅ , … , 1 ⋅ , ⋅ .

The maximum molar concentration is always the same for all places in a SPN while the amount of levels 1 can be different for each place. Hence, both parameters and can be determined globally in the settings component (see Section 6.2) so that they common to all place components and the parameter can also be defined locally in the place components (see Section 6.1.2). For the following explanations, it is assumed that the parameter is defined globally, i.e. it is the same for all places in the SPN.

The derivation of the stochastic level hazard function is shown exemplarily by the following reaction

: 2 3 → ,

where is the deterministic rate constant. The increase of substance can be expressed by mass action kinetics

,

where , , and denote the concentrations of substances , , and . The new concentration can be derived from the current concentration by Euler’s method

∆ ⋅

⇔ ∆ ∆

, with ∆ .

But the abstract concentrations can only increase by one level, i.e. by a concentration of

⁄ , when the reaction occurs, hence,

∆ ∆

,

whereby ∆ is the time that is needed to increase the concentration of substance by one concentration level / . This time is taken as expected value for the occurrence of the reaction and, thus, the hazard function is given by

1

∆ ⋅ ⋅ ,

whereby and are the discrete levels which correspond to the concentrations and and, hence,

⋅ and ⋅ .

The general stochastic level hazard function is then given by

⋅ ⋅

→

∈

, Eq. 5-20

where is the transition-specific deterministic rate constant. The conversion from the deterministic rate constant to the stochastic rate constant can be found in (Wilkinson 2006).

Several methods have been proposed to perform the stochastic simulation of a SPN model either with stochastic mass action hazard functions or stochastic level hazard functions. The first, called Gillespie’s direct method, calculates explicitly which reaction occurs next and

when it occurs (Gillespie 1977). Thereby, reaction occurs with a probability ⁄∑ and the time to the next reaction is the random quantity ∑ . The algorithm is summarized in Algorithm A7 (Appendix A1). This algorithm is direct in the sense that it chooses the reaction and the time when it is occurs directly by generating two random numbers at each iteration.

Gillespie also developed an alternative called first reaction method, which gets along with one random number per iteration. Thereby, a putative time is generated for each reaction.

This is the time when the reaction would occur if no other reaction occured before. Thus, the reaction with the smallest putative time will occur. The method is formalized by Algorithm A8 (Appendix A1).

Gibson and Bruck modified the first reaction method to make it much more effective (Gibson and Bruck 2000). It is called next reaction method. They reduced the number of required random numbers by reusing the . A new putative time is only generated for the reaction that has just occurred. If the hazard function is not changed by this reaction, the remains unaltered. However, if the hazard function is affected by this reaction, the is rescaled appropriately. The knowledge about which hazard function is affected by which reaction is attained from a dependency graph. The method can be performed by the following steps.

Algorithm 1

1. Set the initial numbers of molecules, set 0, calculate the hazard function , and generate a putative time ~ for all .

2. Let the reaction whose putative time is the smallest.

3. Update the number of molecules affected by reaction . 4. Update and generate a new putative time .

5. For each reaction whose hazard function is affected by reaction :

a. Set ,

b. Update ,

c. Set .

6. If stop; otherwise, set and go to step 2.

The next reaction method is used within this study because the dependency graph is already created by the Modelica-tool to perform the hybrid simulation (see Section 3.1.7). In this

manner, the hazard functions are only recalculated when they are changed by firing other transitions. The implementation by the Modelica-language is explained in Section 6.1.3.