ODE to a Petri net

Writing ordinary differential equations (ODEs) to model various natural world processes comes more readily to some than to others.  And once written, it can take some effort to pick apart their real-world meaning.  It would be nice to have a more visual form.  Petri nets are one such approach.

Epidemics

For example, consider the simple SIR model of epidemic infection.

• +++S+++ : those uninfected, but susceptible -- their number is reduced as they become infected, at a rate proportional to the number of susceptibles and the number of infected
• +++I+++ : those infected -- their number is increased by susceptibles who become infected;  it is also reduced as those infected recover, at a rate proportional to the number of infected
• +++R+++ : those recovered -- their number is increased by those who were infected recovering

For those happy with ODEs, it is straightforward to write down a set of coupled equations to model this:

$$\dot{S} = - i S I $$

$$\dot{I} = i S I - r I$$

$$\dot{R} = r I$$
In addition to natural language (the bullet list) and maths (the ODEs), there is another language useful to explain and understand models: diagrams.

For example, we could draw a simple state transition diagram to show the movement from susceptible to infected to recovered:

This captures some of the information, but not all of it.  A (continuous) Petri net can do better:

The circles are called "places", and represent the "things" involved: here the susceptibles, the infected, and the recovered.  The rectangles (the colours aren't significant) are called "transitions", and represent how the things in "input" places get transformed into things in "output" places.

• transition +++r+++: an infected comes in, and a recovered comes out
• transition +++i+++: a susceptible and an infected come in; two infected come out (the two infected outputs are the newly infected, and the original infecter)

This diagram has enough information in it to reproduce the original ODEs.  We have one ODE per place, with the terms being given by the transitions feeding that place.

• place +++R+++. +++R+++ has only one transition feeding it: transition +++r+++.  It is feeding into +++R+++ at a rate proportional to all the inputs to +++r+++: here just +++I+++.  Hence +++\dot{R} \propto I+++.  If we call the constant of proportionality (the rate constant) +++r+++, we get +++\dot{R} = r I+++
• place +++S+++. +++S+++ has only one transition, +++i+++, removing stuff from +++S+++.  It is removing at a rate proportional to all the inputs to +++i+++: here +++S+++ and +++I+++.  Calling the rate constant +++i+++, we get  +++\dot{S} = - i S I+++ (the minus sign is there because we a removing from +++S+++, and it is conventional to keep the rate constants positive)
• place +++I+++.  +++I+++ has two transitions: +++i+++ both feeding it (two input arrows) and removing from it (one output arrow), and +++r+++ removing from it.  We get +++\dot{I} = i S I - r I+++

Thus we have recovered the original equations.

We can write this focussing on the transitions, in a more algorithmic way (algorithms, or pseudo-code, are yet a further language we can use to explain, describe, and define processes).

• initialise the rates of change to each place to be +++0+++.  +++\dot{S}, \dot{I}, \dot{R} := 0+++
• transition +++r+++.  This is removing stuff from place +++I+++ and adding it to +++R+++ at a rate +++rI+++. Update the output place +++R+++ and input place +++I+++ appropriately: +++\dot{I} {-}{=}\ rI+++, +++\dot{R} {+}{=}\ rI+++
• transition +++i+++.  This is removing stuff from place +++I+++ and from place +++S+++ at a rate +++iSI+++ and adding it to +++I+++ at a rate +++2iSI+++ (from the two input arrows).  Hence there is a net input to place +++I+++ at a rate +++iSI+++.  Update the output place +++S+++ and net input place +++I+++ appropriately: +++\dot{I} {+}{=}\ iSI+++, +++\dot{S} {-}{=}\ iSI+++

We can write this as a general algorithm:

  for each place Pi
    Pi_dot := 0
  for each transition Ti
    let Pin = < Pin_1, ... , Pin_n > = list of n places,
                                       one for each input arrow of Ti;
        Pout = list of m places, one for each output arrow of Ti;
        t = Ti x Pin_1 x ... x Pin_n
    for each Pi in Pin
      Pi_dot -= t
    for each Pi in Pout
      Pi_dot += t

So now we have a diagrammatic form, and an ODE form, that are equivalent, and an algorithm to translate on to the other.  This is useful, because we can use them interchangeably, without risk of losing information.  In particular, notice how explanation accompanying the Petri net focusses on what is happening in the transitions, whilst that for the ODE form focusses on what is happening to the places.  Having different forms of explanation can be useful in different circumstances (modelling, communication, modification, validation, calculation, etc).

Catalysis

Although that all works well, the handling of the infecter in the +++i+++ transition seems a bit unnatural: infecter goes in, infecter comes out, resulting in an addition and subtraction of this rate.  The infecter is needed for the transition, and affects the rate of the transition, but is not itself changed by the transition.  In chemistry, this is called a catalyst, and there is some special Petri net syntax for it.  We can draw the SIR Petri net above equivalently as:

Here the dashed arrow means that +++I+++ is a catalyst: it is needed for the transition, but is not consumed by the transition. Hence there is now only one arrow out to +++I+++: since the catalyst wasn't consumed, it doesn't need to be replaced; the remaining single arrow represents the newly infected.

The algorithm needs a little  modification:

  for each place Pi
    Pi_dot := 0
  for each transition Ti
    let Pin = < Pin_1, ... , Pin_n > = list of n places, 
                                       one for each input arrow of Ti;
        Pcat (sublist of Pin) = list of catalytic input arrow of Ti;
        Pout = list of m places, one for each output arrow of Ti
    let t = Ti x Pin_1 x ... x Pin_n
    for each Pi in Pin \ Pcat
      Pi_dot -= t
    for each Pi in Pout
      Pi_dot += t

So the catalytic arrows still contribute to the functional form of the overall rate +++t+++, but not to the changes to the specific places.

Logistic equation

Possibly the simplest bounded growth model in biology is the logistic equation:
$$ \dot{N} = rN(1-N/K)$$where +++r+++ is the growth rate, and +++K+++ is the carrying capacity (so when +++K=N+++, +++\dot{N}=0+++).

This can be drawn as an equivalent Petri net:

• transition +++r+++ (birth): one in, two out
• transition +++rK+++ (competition): two in, one out

These two transitions can also be shown in a simpler catalytic form (if maybe not with the same intuition as before):

• transition +++r+++ (birth): one "catalyses" the birth of the other
• transition +++rK+++ (competition): one "catalyses" the death of the other

Lotka-Volterra predator-prey

The simple predator-prey model, usually cast as rabbits and foxes, has rabbits being born, predated on by foxes to produce more foxes, who then die.  A simplistic version of this might be:

• transition +++b+++ (birth): one rabbit in, two rabbits out
• transition +++d+++ (death): one fox dies
• transition +++p+++ (predation): one fox and one rabbit in, two foxes out

This however has a problem: it has a new fox produced every time a rabbit is eaten.  Real foxes need more food than this to reproduce.  We can't solve the problem by changing the rate +++p+++, as this affects the consumption of rabbits and production of foxes equally.  What we really need is for the consumption of a rabbit to produce a bit of a fox.  We can do this by adding a separate rate to the arrow:

• transition +++p+++ (predation): one fox and one rabbit in, one plus +++\epsilon+++ foxes out

The algorithm needs to be updated to multiply the rate by the weight of the arrow before adding/subtracting  as appropriate.  This then yields the familiar equations:

$$\dot{R} = R(b-pF)$$

$$\dot{F} = F(-d+p\epsilon R)$$

If we use the catalytic form, the diagram simplifies to:

• transition +++p+++ (predation): one fox "catalyses" the transformation of a rabbit into +++\epsilon+++ of a new fox

Lotka-Volterra competition

The simple competition model, usually cast as rabbits and sheep, has rabbits and sheep being born and dying following their own logistic equation, and also competing with each other for resources.

• transition +++rb,rs+++ (birth): one +++x+++ "catalyses" the birth of the next +++x+++
• transition +++rc,sc+++ (competition): one +++x+++ "catalyses" the death of another +++x+++
• transition +++rs+++ (rabbit/sheep competition): one sheep and one rabbit in, a proportion of each out

$$\dot{R} = R(rb-rcS) + (pr-1)rs RS$$

$$\dot{S} = S(sb-scS) + (ps-1)rs RS$$

Combining predator-prey and competition

The competition example has two logistic "subnets", showing how these diagrammatic forms can be readily combined.  So, for example, we could easily add some foxes to the brew:

$$\dot{S} = S(sb-scS) + (ps-1)rs RS$$

$$\dot{R} = R(rb-rcS) + (pr-1)rs RS -pRF$$

$$\dot{F} = F(-d+p\epsilon R)$$

If the foxes also worried the sheep, the diagram would get messier, but would still visually represent the relationships between the different components.

Diagrams v cartoons

Pictures can be very helpful at getting across ideas, but they have their problems if they are ambiguous, incomplete, or otherwise open to misinterpretation.

The (continuous) Petri nets shown here have a formal meaning: they can be translated into equivalent ODEs. They are not informal "cartoons", merely sketching some part of the meaning.  There is an algorithm from diagrams to equations, and it is possible to build a tools allowing the manipulation of diagrams and equations as two different "concrete syntaxes" of the same underlying model.  Which syntax to use depends on what you are doing: it truly is the best of both worlds.

Acknowledgments

• I first came across the formal link between continuous Petri nets and ODEs on Alexi Sharov's web site
• I was reminded of using Petri nets to model population dynamics on reading David Tanzer's guest post on the Azimuth site
• I drew the diagrams in graphviz
• The maths is formatted with LaTeX and displayed using MathJax