Blog - Petri net programming (part 4b)

This page is a blog article in progress, written by David Tanzer. To discuss this article while it’s being written, visit the Azimuth Forum. Please remember that blog articles need HTML, not Markdown.

*guest post by David Tanzer*</i>

In this article, we will describe the “rate equation” for continuous deterministic Petri nets, which gives a formula for the rate at which the the population counts of the species change over time. But before getting into the details, I’d like to outline how we arrive at the notion of continuous deterministic Petri nets, which, on the surface, sound quite different from the stochastic Petri nets that we introduced two articles back.

This model of a continuous, deterministic Petri net is actually derived from the more fundamental model of a discrete, stochastic Petri net, which we introduced in a previous blog article. Think of chemical reactions, for example. To a good approximation, we can view the amounts of the various substances as continuously varying magnitudes (and then view them as real-valued *concentrations* of particles per unit volume), and the processes as continually occurring, at various concentration-dependent rates. But in reality, at a detailed level, there is actually a discrete, integer valued number of each of the particles, and the reaction events take place by separate, discrete, random events.

In a stochastic process, different “runs” of the experiment will produce different sequences of events, and the collection of all possible sequences will be distributed according to a probability distribution. Each sequence of events will lead to a corresponding sequence of states of the system. At every point in time, therefore, there will be a distribution of possible states of the system. We can then take the expected value of the state at each point int time, and form the *sequence of expected states* for the stochastic process.

So, whereas the stochastic process consists of a non-deterministic collection of possible outcomes, all involving discrete integer states (in a stochastic Petri net, discrete population counts), the sequence of expected states is uniquely determined – a deterministic process – involving states with real-valued components – which result from taking the means of the integer-valued components.

Furthermore, under certain statistical assumptions (cite John’s recent paper), as the population counts get large, and the reaction rates increase, a law of large number will kick in, and the variations of the runs of the stochastic process around the mean will decrease (will cluster more tightly around the mean), so that the deterministic sequence of expected states will become a better and better approximation for the observed runs of the stochastic process.

Let’s illustrate with a physical example. Consider an apparatus where we have a winding spiral pipe, with rough and irregular inner walls. At the bottom it opens into a large jar, and at the top we pour a cup containing small pebbles. Each pebble eventually (bounces) makes its way from the top of the pipe to the bottom and into the jar, but due to the irregularities of friction, etc., there may be a substantial variation in the time that each pebble takes to traverse the pipe.

The events are the arrivals of the pebbles in the jar, and the state at each point we will take to be the number of pebbles in the jar. If we graph the state (the pebble count) as a function of time, it will be a non-decreasing function, with integer values, and integer jumps as pebbles enter the jar. But the *expected* number of pebbles in the jar, as a function of time, will be an increasing, continuous function.

(We suppose the pipe is wide in comparison to the pebble sizes, so there are no clogging or queueing effects.)

Now consider what happens if, say, we multiply by ten the number of pebbles that we release. The results of the process will still be stochastic, but the overall process, as recorded by the rising count of pebbles in the jar, will be smoother, and closer to the averaged out sequence of expected pebble counts.

In fact, by increasing the pebble counts within each run of the experiment, we are performing a kind of averaging there: the effects of the slower pebbles will be offset by the effects of the faster pebbles – the count in the jar reflects the aggregate effects of all the pebbles.

In the limiting process of a large number of pebbles, the behavior looks more like pouring sand into the pipe, which begins to look like a continuous fluid. The continuous fluid model is an even better approximation when we reach the scale of chemical reactions, which involve moles of chemical entities. As we pass to the theoretical limit, involving an infinite number of infinitesimal particles, we reach the fluid model for a continuous, deterministic process.

Now we are working in the continuous flow model for a reaction network.

Let’s first look at the behavior for a reaction network that has just one transition. The transition is like a little pump, that drains fluid from its input containers, and emits fluid on its output containers. The rate at which the pump is running is proportional to the product of the amounts of fluid at each of its input containers. This is the mass action kinetics, which we described two articles back. The explanation for this law, in terms of a fluid of particles, is that a given amount of substance at one of the inputs to a reaction represents a given number of candidate particles for the reaction, and hence will contribute linearly to the probability of a reaction occurring within a small time interval, and hence will contribute linearly to the firing rate of the reaction.

Let’s put a formula to it.

Let Ci be the ith container in the network. Let ai(t) be the amount of substance stored at the ith container, at time t.

Let R be one of the reactions in the network.

As we said before, the firing rate of R at time t is proportional to the amounts of each of its inputs ((times the rate constant for R)):

firingRate(R,t) = rateConstant(R) * Prod {ai(t) | i in INPS(R)}.

Let INPS(R) = the multiset of indicies i in {1,…,n} for the input connections to R.

(Note that if Ci is input multiple times to R, then i will appear multiple times in INPS(R), and the factor ai(t) will appear multiple times in this product.)

Next, what is the effect of the firing rate of R, at time t, on the amount ai(t) stored a container i at that time?

Well, each firing event of R causes one unit of fluid to be drained from each of the input containers, and one unit of fluid to be added to each of the output containers. So, when R has one output to Ci, then ((the rising fluid at Ci due to R is expressed by)):

d/dt ai(t) = firingRate(R,t)

And when R has in input to Ci, then ((The draining fluid at Ci due to R is expressed by)):

d/dt ai(t) = -firingRate(R,t).

In the general case, R may have multiple inputs to Ci, and also, possibly at the same time, multiple outputs to Ci – so the effect of these drainings and infusions must be netted out.

The net change at Ci, per firing of R, is then:

OUT(R)(i) - INP(R)(i).

This gives us:

dai(t) / dt = FiringRate(R,t) (OUT(R)(i) - INP(R)(i)) =

RateConstant(R) * (PROD j in INP(R) aj(t)) * (OUT(R)(i) - INP(R)(i))

This is the rate equation, for a single transition.

The extension to the case of multiple transitions is straightforward, because at a given instant, the firings of each of the transitions are separately determined, as per the mass action law that we just described, and the contributions of each transition to the increasing or decreasing of the fluid at a given container, are addititve – we have have to sum the rates of change given by the equation for a single transition:

dai(t) / dt = SUM(R in REACTIONS) (PROD j in INP(R) aj(t)) * (OUT(R)(i) - INP(R)(i)).

This is the general rate equation.

This is a system of differential equations, one for each container i in the network.

Let’s write a(t) = (a1(t),…,an(t)) for the vector of all amounts stored in the containers, i.e., the complete state of the network.

((Explain how INP(R) can be both a multiset, and a function taking containers to natural numbers.))

We can then write the whole system of equations as a single equation involving vectors:

da(t) / dt = SUM (R in REACTIONS) RateConstant(R) * (PROD j in INP(R) aj(t)) * (OUT(R) - INP(R))

Here, OUT(R) and INP(R) designate the vectors:

OUT(R) = (OUT(R)(1), …, OUT(R)(n)) INP(R) = (INP(R)(1), …, INP(R)(n)),

which indicate the number of times that R outputs or inputs to each of the containers.

The difference OUT(R) - INP(R) is the vector that gives the net increase at each container, due to a single firing of the reaction R.

Looking at this equation, we see that it gives the direction and rate of change of the state of the system, by the vector of derivatives

da(t) / dt = (da1(t) / dt, …, dan(t) / dt),

which is put as a function of the current state of the system a(t) = (a1(t), …, an(t)); all the other terms on the right hand side of the equation, i.e., RateConstant(R), INP(R), and OUT(R), are constants that are fixed for a given reaction network. In a word, this equation gives the derivative of the system as a function of the state of a system.

Therefore, starting out at an initial state a(0), this equation will determine a unique trajectory for the system, as a function of time. (This is a theorem that follows from the structure of the differential equation.) It then gives a deterministic law of motion for the system – this is the general meaning of the rate equation.

(And solving for da(t) / dt = 0 gives the conditions for system equilibrium.)

in the next article, we will look at specific examples of reaction networks, and apply the rate equation to them, to determine the structure of the evolution of the system, and their states of equilibrium. Then, for simple but basic systems, we will give graphical “streamlines” that show the evolution of the systems.

Let’s consider an important, basic example: a pair of opposite reactions.

r: a1 –alpha–> a2 s: a2 –beta–> a1

INP(r) = {1} = (1,0) OUT(r) = {2} = (0,1) INP(s) = {2} = (0,1) OUT(s) = {1} = (1,0)

OUT(r) - INP(r) = (-1,1) OUT(s) - INP(s) = (1,-1)

RateConstant(r) = alpha RateConstant(s) = beta

da(t) / dt = (da1(t) / dt, da2(t) / dt)

RateConstant(r) * PROD(j in {1}) aj(t)) (-1,1) +

RateConstant(s) * PROD(j in {2}) aj(t)) (1,-1) =

alpha a1(t) (-1, 1) + beta a2(t) (1, -1) =

(- alpha a1(t) + beta a2(t), alpha a1(t) - beta a2(t))

Note that da1 / dt = - da2 / dt,

which makes sense, because whatever tokens are being lost on the one container must be being gained on the other.

Equilibrium is reached when da / dt = 0, which means that

alpha a1(t) = beta a2(t),

so a1/a2 = beta / alpha. Thus, beta/alpha gives the proportions of a1 and a2 at equilibrium.

category: blog