Petri.jl version of SIR model demonstrating conversion to ODE

ode_petri.jmd

# Ordinary differential equation model using Petri.jl
Simon Frost (@sdwfrost), 2020-07-08

## Introduction


This implementation considers the SIR model as a Petri net, based on the example of the [documentation of `Petri.jl`](https://mehalter.github.io/Petri.jl/stable/usage/), which is then used to generate an ordinary differential equation model.

## Libraries

```julia
using DifferentialEquations
using OrdinaryDiffEq
using Petri
using LabelledArrays
using DataFrames
using Plots
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph
```

## Transitions

The Petri model is specified using a vector of the model states (as symbols), and a labelled vector of the transition rates; in this case, `inf` (infection) and `rec` (recovery). Each transition is a tuple of labeled vectors with inputs and outputs.

```julia
sir = Petri.Model([:S,:I,:R],LVector(
                                inf=(LVector(S=1,I=1), LVector(I=2)),
                                rec=(LVector(I=1),     LVector(R=1))))
```

The resulting Petri net can be converted into an ordinary differential equation model as follows.

```julia
sir_ode! = vectorfields(sir)
```

Using Graphviz, a graph showing the states and transitions can also be generated from the Petri net.

```julia
graph = Graph(sir)
```

## Time domain

We set the timespan for simulations, `tspan`.

```julia
tspan = (0.0,40.0)
```

## Initial conditions

As the model is given in terms of labelled vectors, we have to specify the initial conditions the same way.

```julia
u0 = LVector(S=990.0,I=10.0,R=0.0)
```

## Parameter values

Define the parameters of the model, `p`. For the Petri net model, each parameter corresponds to the rate of a transition.

```julia
p0 = LVector(β=0.05,c=10.0,γ=0.25,N=1000.0)
p = LVector(inf=p0.β*p0.c/p0.N,rec=p0.γ)
```

## Running the model

```julia
prob_ode = ODEProblem(sir_ode!,u0,tspan,p)
```

As `Petri.jl` has its own `solve` method, we have to specify that we want to run `solve` from `OrdinaryDiffEq`.

```julia
sol_ode = OrdinaryDiffEq.solve(prob_ode,Tsit5());
```

## Plotting

We can now plot the results.

```julia
plot(sol_ode)
```

```{julia; echo=false; skip="notebook"}
include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()
```