*guest post by David A. Tanzer*

### An introduction to stochastic Petri nets

In the previous article, I explored a simple computational model called Petri nets. They are used to model reaction networks, and have applications in a wide variety of fields, including population ecology, gene regulatory networks, and chemical reaction networks. I presented a simulator program for Petri nets, but it had an important limitation: the model and the simulator contain no notion of the *rates* of the reactions. But these rates critically determine the character of the dynamics of network.

Here I will introduce the topic of ‘stochastic Petri nets,’ which extends the basic model to include reaction dynamics. Stochastic means random, and it is presumed that there is an underlying *random process* that drives the reaction events. This topic is rich in both its mathematical foundations and its practical applications. A direct application of the theory yields the rate equation for chemical reactions, which is a cornerstone of chemical reaction theory. The theory also gives algorithms for analyzing and simulating Petri nets.

We are now entering the ‘business’ of software development for applications to science. The business logic here is nothing but math and science itself. Our study of this logic is not an academic exercise that is tangential to the implementation effort. Rather, it is the first phase of a *complete* software development process for scientific programming applications.

The end goals of this series are to develop working code to analyze and simulate Petri nets, and to apply these tools to informative case studies. But we have some work to do *en route*, because we need to truly understand the models in order to properly interpret the algorithms. The key questions here are when, why, and to what extent the algorithms give results that are empirically predictive. We will therefore be embarking on some exploratory adventures into the relevant theoretical foundations.

The overarching subject area to which stochastic Petri nets belong has been described as *stochastic mechanics* in the network theory series here on Azimuth. The theme development here will partly parallel that of the network theory series, but with a different focus, since I am addressing a computationally oriented reader. For an excellent text on the foundations and applications of stochastic mechanics, see:

• Darren Wilkinson, *Stochastic Modelling for Systems Biology*, Chapman and Hall/CRC Press, Boca Raton, Florida, 2011.

### Review of basic Petri nets

A Petri net is a graph with two kinds of nodes: species and transitions. The net is populated with a collection of ‘tokens’ that represent individual entities. Each token is attached to one of the species nodes, and this attachment indicates the type of the token. We may therefore view a species node as a container that holds all of the tokens of a given type.

The transitions represent conversion reactions between the tokens. Each transition is ‘wired’ to a collection of input species-containers, and to a collection of output containers. When it ‘fires’, it removes one token from each input container, and deposits one token to each output container.

Here is the example we gave, for a simplistic model of the formation and dissociation of H_{2}O molecules:

The circles are for species, and the boxes are for transitions.

The transition **combine** takes in two H tokens and one O token, and outputs one H_{2}O token. The reverse transition is **split**, which takes in one H_{2}O, and outputs two H’s and one O.

An important application of Petri nets is to the modeling of *biochemical* reaction networks, which include the gene regulatory networks. Since genes and enzymes are molecules, and their binding interactions are chemical reactions, the Petri net model is directly applicable. For example, consider a transition that inputs one gene G, one enzyme E, and outputs the molecular form G • E in which E is bound to a particular site on G.

Applications of Petri nets may differ widely in terms of the population sizes involved in the model. In general chemistry reactions, the populations are measured in units of moles (where a mole is ‘Avogadro’s number’ 6.022 · 10^{23} entities). In gene regulatory networks, on the other hand, there may only be a handful of genes and enzymes involved in a reaction.

This difference in scale leads to a qualitative difference in the modelling. With small population sizes, the stochastic effects will predominate, but with large populations, a continuous, deterministic, average-based approximation can be used.

### Representing Petri nets by reaction formulas

Petri nets can also be represented by formulas used for chemical reaction networks. Here is the formula for the Petri net shown above:

H_{2}O ↔ H + H + O

or the more compact:

H_{2}O ↔ 2 H + O

The double arrow is a compact designation for *two* separate reactions, which happen to be opposites of each other.

By the way, this reaction is not physically realistic, because one doesn’t find isolated H and O atoms traveling around and meeting up to form water molecules. This is the actual reaction pair that predominates in water:

2 H_{2}O ↔ OH^{-} + H_{3}O^{+}

Here, a hydrogen nucleus H^{+}, with one unit of positive charge, gets removed from one of the H_{2}O molecules, leaving behind the **hydroxide** ion OH^{-}. In the same stroke, this H^{+} gets re-attached to the other H_{2}O molecule, which thereby becomes a **hydronium** ion, H_{3}O^{+}.

For a more detailed example, consider this reaction chain, which is of concern to the ocean environment:

CO_{2} + H_{2}O ↔ H_{2}CO_{3} ↔ H^{+} + HCO_{3}^{-}

This shows the formation of **carbonic acid**, namely H_{2}CO_{3}, from water and carbon dioxide. The next reaction represents the splitting of carbonic acid into a hydrogen ion and a negatively charged **bicarbonate ion**, HCO_{3}^{-}. There is a further reaction, in which a bicarbonate ion further ionizes into an H^{+} and a doubly negative **carbonate ion** CO_{3}^{2-}. As the diagram indicates, for each of these reactions, a reverse reaction is also present. For a more detailed description of this reaction network, see:

• Stephen E. Bialkowski, Carbon dioxide and carbonic acid.

Increased levels of CO_{2} in the atmosphere will change the balance of these reactions, leading to a higher concentration of hydrogen ions in the water, i.e., a more acidic ocean. This is of concern because the metabolic processes of aquatic organisms is sensitive to the pH level of the water. The ultimate concern is that entire food chains could be disrupted, if some of the organisms cannot survive in a higher pH environment. See the Wikipedia page on ocean acidification for more information.

**Exercise.** Draw Petri net diagrams for these reaction networks.

### Motivation for the study of Petri net dynamics

The relative rates of the various reactions in a network critically determine the qualitative dynamics of the network as a whole. This is because the reactions are ‘competing’ with each other, and so their relative rates determine the direction in which the state of the system is changing. For instance, if molecules are breaking down faster then they are being formed, then the system is moving towards full dissociation. When the rates are equal, the processes balance out, and the system is in an equilibrium state. Then, there are only temporary fluctuations around the equilibrium conditions.

The rate of the reactions will depend on the number of tokens present in the system. For example, if any of the input tokens are zero, then the transition can’t fire, and so its rate must be zero. More generally, when there are few input tokens available, there will be fewer reaction events, and so the firing rates will be lower.

Given a specification for the rates in a reaction network, we can then pose the following kinds of questions about its dynamics:

• Does the network have an equilibrium state?

• If so, what are the concentrations of the species at equilibrium?

• How quickly does it approach the equilibrium?

• At the equilibrium state, there will still be temporary fluctuations around the equilibrium concentrations. What are the variances of these fluctuations?

• Are there modes in which the network will oscillate between states?

This is the grail we seek.

Aside from actually performing empirical experiments, such questions can be addressed either analytically or through simulation methods. In either case, our first step is to define a theoretical model for the dynamics of a Petri net.

### Stochastic Petri nets

A **stochastic Petri net** (with kinetics) is a Petri net that is augmented with a specification for the reaction dynamics. It is defined by the following:

• An **underlying Petri net**, which consists of species, transitions, an input map, and an output map. These maps assign to each transition a multiset of species. (Multiset means that duplicates are allowed.) Recall that the state of the net is defined by a *marking* function, that maps each species to its population count.

• A **rate constant** that is associated with each transition.

• A **kinetic model**, that gives the expected firing rate for each transition as a function of the current marking. Normally, this kinetic function will include the rate constant as a multiplicative factor.

A further ‘sanity constraint’ can be put on the kinetic function for a transition: it should give a positive value if and only if all of its inputs are positive.

• A **stochastic model**, which defines the probability distribution of the time intervals between firing events. This specific distribution of the firing intervals for a transition will be a function of the expected firing rate in the current marking.

This definition is based on the standard treatments found, for example in:

• M. Ajmone Marsan, Stochastic Petri nets: an elementary introduction, in *Advances in Petri Nets*, Springer, Berlin, 1989, 1–23.

or Wilkinson’s book mentioned above. I have also added an explicit mention of the kinetic model, based on the ‘kinetics’ described in here:

• Martin Feinberg, Lectures on chemical reaction networks.

There is an implied *random process* that drives the reaction events. A classical random process is given by a container with ‘particles’ that are randomly traveling around, bouncing off the walls, and colliding with each other. This is the general idea behind Brownian motion. It is called a random process because the outcome results from an ‘experiment’ that is not fully determined by the input specification. In this experiment, you pour in the ingredients (particles of different types), set the temperature (the distributions of the velocities), give it a stir, and then see what happens. The outcome consists of the paths taken by each of the particles.

In an important limiting case, the stochastic behavior becomes deterministic, and the population sizes become continuous. To see this, consider a graph of population sizes over time. With larger population sizes, the relative jumps caused by the firing of individual transitions become smaller, and graphs look more like continuous curves. In the limit, we obtain an approximation for high population counts, in which the graphs are continuous curves, and the concentrations are treated as continuous magnitudes. In a similar way, a pitcher of sugar can be approximately viewed as a continuous fluid.

This simplification permits the application of continuous mathematics to study of reaction network processes. It leads to the basic *rate equation* for reaction networks, which specifies the direction of change of the system as a function of the current state of the system.

In this article we will be exploring this continuous deterministic formulation of Petri nets, under what is known as the mass action kinetics. This kinetics is one implementation of the general specification of a kinetic model, as defined above. This means that it will define the *expected* firing rate of each transition, in a given marking of the net. The probabilistic variations in the spacing of the reactions—around the mean given by the expected firing rate—is part of the stochastic dynamics, and will be addressed in a subsequent article.

### The mass-action kinetics

Under the **mass action kinetics**, the expected firing rate of a transition is proportional to the product of the concentrations of its input species. For instance, if the reaction were A + C → D, then the firing rate would be proportional to the concentration of A times the concentration of C, and if the reaction were A + A → D, it would be proportional to the square of the concentration of A.

This principle is explained by Feinberg as follows:

For the reaction A+C → D, an occurrence requires that a molecule of A meet a molecule of C in the reaction, and we take the probability of such an encounter to be proportional to the product [of the concentrations of A and C]. Although we do not presume that every such encounter yields a molecule of D, we nevertheless take the occurrence rate of A+C → D to be governed by [the product of the concentrations].

For an in-depth proof of the mass action law, see this article:

• Daniel Gillespie, A rigorous definition of the chemical master equation, 1992.

Note that we can easily pass back and forth between speaking of the population counts for the species, and the concentrations of the species, which is just the population count divided by the total volume V of the system. The mass action law applies to both cases, the only difference being that the constant factors of (1/V) used for concentrations will get absorbed into the rate constants.

The mass action kinetics is a basic law of empirical chemistry. But there are limits to its validity. First, as indicated in the proof in the Gillespie, the mass action law rests on the assumptions that the system is well-stirred and in thermal equilibrium. Further limits are discussed here:

• Georg Job and Regina Ruffler, *Physical Chemistry* (first five chapters), Section 5.2, 2010.

They write:

…precise measurements show that the relation above is not strictly adhered to. At higher concentrations, values depart quite noticeably from this relation. If we gradually move to lower concentrations, the differences become smaller. The equation here expresses a so-called “limiting law“ which strictly applies only when c → 0.

In practice, this relation serves as a useful approximation up to rather high concentrations. In the case of electrically neutral substances, deviations are only noticeable above 100 mol m

^{−3}. For ions, deviations become observable above 1 mol m^{−3}, but they are so small that they are easily neglected if accuracy is not of prime concern.

Why would the mass action kinetics break down at high concentrations? According to the book quoted, it is due to “molecular and ionic interactions.” I haven’t yet found a more detailed explanation, but here is my supposition about what is meant by molecular interactions in this context. Doubling the number of A molecules doubles the number of expected collisions between A and C molecules, but it also reduces the probability that any given A and C molecules that are within reacting distance will actually react. The reaction probability is reduced because the A molecules are ‘competing’ for reactions with the C molecules. With more A molecules, it becomes more likely that a C molecule will simultaneously be within reacting distance of several A molecules; each of these A molecules reduces the probability that the other A molecules will react with the C molecule. This is most pronounced when the concentrations in a gas get high enough that the molecules start to pack together to form a liquid.

### The equilibrium relation for a pair of opposite reactions

Suppose we have two opposite reactions:

Since the reactions have exactly opposite effects on the population sizes, in order for the population sizes to be in a stable equilibrium, the expected firing rates of and must be equal:

By mass action kinetics:

where means the concentration of

Hence at equilibrium:

So:

where is the equilibrium constant for the reaction pair.

### Equilibrium solution for the formation and dissociation of a diatomic molecule

Let A be some type of atom, and let D = A_{2} be the diatomic form of A. Then consider the opposite reactions:

From the preceding analysis, at equilibrium the following relation holds:

Let and be the population counts for A and B, and let

be the total number of units of A in the system, whether they be in the form of atoms or diatoms.

The value of is an invariant property of the system. The reactions cannot change it, because they are just shuffling the units of A from one arrangement to the other. By way of contrast, is not an invariant quantity.

Dividing this equation by the total volume , we get:

where is the concentration of the units of A.

Given a fixed value for and the rate constants and , we can then solve for the concentrations at equilibrium:

Since can’t be negative, only the positive square root is valid.

Here is the solution for the case where :

### Conclusion

We’ve covered a lot of ground, starting with the introduction of the stochastic Petri net model, followed by a general discussion of reaction network dynamics, the mass action laws, and calculating equilibrium solutions for simple reaction networks.

We still have a number of topics to cover on our journey into the foundations, before being able to write *informed* programs to solve problems with stochastic Petri nets. Upcoming topics are (1) the deterministic rate equation for general reaction networks and its application to finding equilibrium solutions, and (2) an exploration of the *stochastic* dynamics of a Petri net. These are the themes that will support our upcoming software development.

Great post! For info, the correct reference for the first edition of my book is:

Wilkinson, D. J. (2006) Stochastic modelling for systems biology, Boca Raton, Florida: Chapman and Hall/CRC Press.

but people may prefer the second edition:

Wilkinson, D. J. (2011) Stochastic modelling for systems biology, second edition, Boca Raton, Florida: Chapman and Hall/CRC Press.

especially since this edition considers in more detail the problem of developing software for simulation of stochastic Petri nets, and describes a free open source software package for it. Further information can be found at the web site:

http://www.staff.ncl.ac.uk/d.j.wilkinson/smfsb/2e/

David wrote:

Hey, you just made me realize: we can use stochastic Petri nets obeying the law of mass action to study the limits of this law’s validity! For example, say we’ve got a reaction

A + B → C

Then we can make up a stochastic Petri net describing a system with n boxes, and thus 3n types of particle: A

_{i}, B_{i}, and C_{i}for 1 ≤ i ≤ n. We include reactionsA

_{i}+ B_{i}→ C_{i}all with the same rate constant, along with reactions that allow particles to hop from a box to any neighboring box. We can define ‘neighboring box’ in various ways, and choose whatever rate constants we like for these ‘hopping’ reactions.

We should see some interesting effects, since each individual box will count as ‘well-stirred’, but the whole system will not be well-stirred unless the rate constants for hopping from one box to another are high.

That sounds like an island model or stepping stone model in population genetics. Imagine a population of birds living on an island which is then separated by sea level rise into two islands. Mostly the birds stay on the island that they were born on and exchange alleles with one another, but occasionally they migrate to the other. How much migration between the two islands is necessary to prevent the two populatons drifting apart (genetically speaking) and eventually becoming two species? A rule of thumb is that one migration per generation is enough.

Here’s an old reference from 1971

http://www.genetics.org/content/73/1/147

and a recent review

http://dna.ac/filogeografia/PDFs/popgen/Charlesworth_09_Effective_Population_Size.pdf

Good point! By the way, at that biodiversity conference I went to, Lou Jost said (quite convincingly) that the ‘one migration per generation’ rule is completely wrong—an artifact of very silly ways of measuring biodiversity! The rule seems counterintuitive, and there’s a reason why: it’s wrong.

I’d guess that the “silly ways of measuring biodiversity” are using entropic measures (e.g., the Shannon index) instead of effective-number ones. Was Jost’s critique something along those lines?

Blake wrote:

Yes. I forget precisely which biodiversity measure was used to derive this ‘one migration per generation’ rule, but it wasn’t an effective number, and I think Jost or someone showed that if you use an effective number this rule changes to something more plausible.

(It’s not plausible that for two arbitrarily large populations of birds on separate islands, one bird flying from one island to the other per generation is enough to keep them from genetically diverging and becoming separate species!)

As a Post Keynesian economist, I am interested in how weighted Petri nets could be used to simulate processes of economic production and reproduction along Classical rather than Neoclassical lines. Kurz and Salvadori draw a useful distinction between neoclassical modelling and the multi-sectoral modelling approach developed by Sraffa, Leontieff, and Von Neumann. The latter two researchers also made notable contributions to the theory of automatons. To represent Classical systems, Von Neumann’s input-output systems (which can readily be represented by Petri nets) would have to be augmented with cost information for each transition (along with time-related storage costs for tokens at each place). These costs would have to include rates of profit for each industry. A secondary process would also have to be modeled – the equalization of the rate of profit across each sector achieved through the flow of capital out of low profit sectors into high profit sectors. Sraffa argued that Marx was constrained through lack of access to linear algebra. You could also claim that Sraffa was constrained by lack of access to process algebra or coalgebra! Your thoughts on this would be much appreciated.

Hi! Your idea sounds interesting, but I don’t know enough economics to understand in detail what you’re talking about. Can you describe a tiny example of the kind of model you’re imagining?

I’d also enjoy seeing an example of one of von Neumann’s ‘input-output systems’. Presumably there are examples lurking around the web.

Hi John! References to Leontief & von Neumann’s respective input-output modeling approaches can be found in G. Bonnano’s 1995 paper “Modelling production with Petri Nets” on his website: http://www.econ.ucdavis.edu/faculty/bonanno/PDF/petri.pdf

And references to the same authors in regard to automatons can be found in P. Mirowski & K. Somefun’s 1998 paper “Towards an Automata approach of Institutional (and Evolutionary) Economics” available at:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.131.7417&rep=rep1&type=pdf

But no one I know of has modelled Sraffa’s (1960) “Production of Commodities by Means of Commodities” using coalgebras. Linear input-output models were the stock in trade of development modeling in the 1950 but have been largly “superceded” by Computerised General Equilibrium models. Many National Statistics Bureaus around the world publish I-O accounts.

Thanks for the pointer to Bonnano’s paper. It seems he only uses Petri nets to say how certain outputs can be made from certain inputs, not to model

dynamics, that is, how things change with time. Of course for dynamic one wants something more like a stochastic Petri net. Without dynamics one tends to focus on questions akin to the ‘reachability problem’: that is, starting with some collection of things, which other things can you eventually make? There’s been a huge amount of work on this in the Petri net literature, which I briefly sketched here.With my background in physics, I’m drawn toward questions of dynamics.

I think I’ve heard of Leontief models being used to make economic predictions… but quickly reading a bit abou tthis, I get the feeling there’s no true dynamics here, just a linear map sending a vector of inputs to a vector of outputs, which could let you predict the amount of outputs if you knew the vector of inputs (if the model were right). Is that right?

There must be a huge hunger for going beyond ‘equilibrium’ models, no? I’m imagining a world where weather forecasters only had equilibrium models, and thinking about how unsatisfactory this would be.

I don’t know his work, either. Could you say a word about why coalgebras could be relevant?

Wikipedia has a good entry on Sraffa. In response to your thoughful comments, let me first note that Wassily Leontief did develop dynamic versions of his input-output models.

Second, I would like to point to the fact that Piero Sraffa included in the title of his book “Prelude to a Critique of Economic Theory”. This work, together with earlier complaints by Joan Robinson, launched the debates on capital theory. Amongst other things, Sraffa shows that in a multisectoral economy, with differing ratios of dead to living labour in each sector, there will be no simple monotonic functional relationship between the rate of profit and the capital intensity of chosen production techniques.

Sraffa used linear algebra and the Perron-Frobenius theorems to establish his findings, but there are different orders of temporality at play in his model: capital flows between sectors to equalize the rate of profit while economies and diseconomies of scale come into play. These processes are clearly longer-term, while changes in the distribution of income between wages and profits are of a shorter-term nature. This is exactly where a coalgebraic approach could come into play.

Here, I am not thinking so much of hybrid systems (i.e. discrete versus continuous) but more a case where you have a nesting of automatons or Petri nets within other nets to accommodate processes operating in accordance with different orders of temporality. Needless to say, this will complicate things significantly in regard to standard measures of liveliness, reachability, safety etc.

As Mirowski points out, most applications of automatons in economics have been designed to account for such things as bounded rationality in a game theory context, or in his case, for the evolution of different market forms: single-auction, double-auction etc! Instead, I want to use process algebraic techniques to remodel dynamics processes at work in the Classical economics of Ricardo and Marx.

James, I found your post to be interesting, and it was good to read the Bonanno paper — which opens the window to another application area of Petri nets.

I don’t have enough knowledge of economics to make a substantive assessment of the feasibility of your idea.

That given, I did make an attempt to imagine the kind of Petri net you are describing. I didn’t get very far, but here is the thought process:

What are the places in this economic Petri net? In the Bonanno paper, production processes are described by transitions, and the places correspond to definite commodities in the network — each transition consumes a certain amount of its input commodities, and generates a certain amount of its output commodities.

Now you want to add the profit generated by the firing of each transition. So in a first approximation, that would be an additional output arc, that adds a quantity of profit into some accumulation node. These nodes would be…bank accounts of the owners of the process? That would mean using many micro-instances of a transition, say one for each factory owner. Or could it be an aggregate profit account.

But this raises the questions of modeling circulation. The profit will get reinvested in more commodities, and how would that be represented in the Petri net.

Next, you want to introduce sectoral rates of profit. Will these be places in the Petri net? Each firing of a transition will then result in a variable amount of tokens being transferred to the output account — the product of the sector rate times the input value. How is input value to be represented in the net? Are there places representing prices in the net. But then how do we accomplish multiplications in a Petri net, e.g., the value of input A is the product of the number of tokens there times the price of A. Similarly for the multiplication by the sectoral rate of profit.

This suggests that a more powerful model of computation may be needed, closer to a general dataflow graph — say something like where the firing of a transition leads to the evaluation of a formula, whose results then get deposited at the output places.

Next, I haven’t even begun to picture how the equalization of sectoral profits would be represented by the net. Can you give an references to a good computational model for this equalization process — that could serve as a point of departure for thinking about Petri-netification of this part of the process.

Anyone care to make another foray into this topic, maybe you can get further than I did here, or add some clarification.

Or do we simplify, and use monetary units for all of the commodities…

Hi David, and thanks for a thoughtful deliberation on what I am proposing. I have read some literature on related matters such as Adulla & Mayr on “Minimal Cost Reachability/Coverability in Priced

Timed Petri Nets” who propose using priced timed Petri nets, and Gradisar & Music on “Automated Petri-net modelling based on production management data”, who construct a timed Petri net simulator using their own MatLab-based software routines.

However, as you point out, long-run equalization of the rate of profit would have to be achieved using some kind of stock adjustment (difference equation) process rather than a Petri net so we are looking at a kind of hybrid system (though one that is still discrete and recursive). You can see why I have been grappling with this modelling issue for some time.

Gradual adjustment also means that each sector would be using a different rate of profit for the spot pricing of its outputs and this would feed in to other users of these outputs(with each represented by individual transitions) as inputs into their own production.

Finally, the presence of economies or diseconomies of scale would also influence unit costs of production! This also means, as Sraffians recognize, that the Perron-Frobenius conditions on the input-output matrix will change with shifts in the composition of output that are induced by changes in income distribution. Even a linear multi-sectoral economy is a pretty complex entity!

I suppose this would also upset John’s ‘quantum stochastics’ solution to stochatic Petri nets (which use Matrix exponentiation and the Perron-Frobenius condition). Any further suggestions on your part would be much appreciated!

[...] Last time we looked at stochastic Petri nets, which use a random event model for the reactions. Individual entities are represented by tokens that flow through the network. When the token counts get large, we observed that they can be approximated by continuous quantities, which opens the door to the application of continuous mathematics to the analysis of network dynamics. [...]