diff --git a/docs/literate/covid/disease_strains.jl b/docs/literate/covid/disease_strains.jl index a9592217..749c1f37 100644 --- a/docs/literate/covid/disease_strains.jl +++ b/docs/literate/covid/disease_strains.jl @@ -10,10 +10,11 @@ using DisplayAs, Markdown # This example presents models incorporating multiple strains of disease and vaccine type. # Importantly, it shows why stratification by disease strain is different from other stratifications, e.g. geography or age, and requires using a different type system. +# If you are unfamiliar with stratification, we reccomend first reading the [stratification tutorial](@ref stratification_example). # ## Define basic epidemiology model -# We start by defining our basic type system for infectious disease models. +# We start by defining our basic type system, $P_{infectious}$ for infectious disease models. const infectious_ontology = LabelledPetriNet( [:Pop], @@ -24,9 +25,9 @@ const infectious_ontology = LabelledPetriNet( to_graphviz(infectious_ontology) -# We define a simple SIRD model with reflexive transitions typed as `:strata` to indicate which states can be stratified. -# Here we add reflexive transitions to the susceptible, infected, and recovered populations but we leave out the dead -# population because they cannot do things such as get vaccinated or travel between regions. +# We define a simple SIRD model with reflexive transitions typed as "strata" to indicate which stratified +# states will interact with transitions amongst strata. All but the dead population will have reflexive +# transitions amongst strata because no individuals may leave any death state. sird_uwd = @relation (S,I,R,D) where (S::Pop, I::Pop, R::Pop, D::Pop) begin infect(S, I, I, I) @@ -41,9 +42,10 @@ to_graphviz(dom(sird_model)) # ## Define a model of multiple vaccine types -# We also define a model of vaccination with multiple vaccine types. -# In this model, vaccination transitions are typed as `:strata`. -# Note that the `:infect` transitions must be included to enable cross-infection between different vax types. +# We also define a model of vaccination with multiple vaccine types, typed over $P_{infectious}$. +# Each stratum has reflexive "disease" transitions. Transitions which represent the administration +# of a vaccine to an individual are typed as "strata". Infection (typed as "infect") between individuals of different +# vaccination status is possible (i.e.; we do not assume a perfect vaccine). function vax_model(n) uwd = RelationDiagram(repeat([:Pop], n+1)) @@ -129,9 +131,10 @@ to_graphviz(dom(strain_model′(2))) # ## Stratify the SIRD model for two strains -# Unfortunately, stratification of these models does not produce the desired result. -# There are quite a few extraneous states and transitions. -# The primary issue is the asymmetry in the role of the uninfected population. +# Unfortunately, stratification of these models does not produce the desired result, leading +# to many states and transitions. The problem is that it does not make sense for the +# uninfected population to be stratified over the full set of states of the SIRD model +# (i.e.; uninfected persons can never have any disease state other than "S"). # We can address this by changing the type system. typed_product(sird_model, strain_model′(2)) |> dom |> to_graphviz @@ -153,7 +156,8 @@ const strain_ontology = LabelledPetriNet( to_graphviz(strain_ontology) -# We now reform the SIRD model using the new type system. +# We now reform the SIRD model using the new type system. Note that the `where` clause is used to define +# the types for boxes in the UWD. sird_for_strains_uwd = @relation (S,I,R,D) where (S::Uninf, I::Inf, R::Inf, D::Inf) begin infect(S, I, I, I) disease(I, R) @@ -194,23 +198,29 @@ end to_graphviz(dom(strain_model(2))) -# When we now stratify we get the desired model. +# When we now stratify we get the desired model. Note however that there are extraneous "strata" +# transitions; this is because there were no non-trivial stratum transitions in either $\phi$ or $\phi^{'}$ +# (note that only "disease" and "infect" transition types were needed to fully define the multiple strain dynamics). sird_strain = typed_product(sird_for_strains_model, strain_model(2)) to_graphviz(dom(sird_strain)) -# ## Post-composition: Typing the type system +# ## Post-composition: Retyping the type system # In some instances, we may want to relate models typed to different type systems. # For example, we usually type our `simple_trip` model of geographic regions to the `infectious_ontology` such that we can stratify a disease model by geographic regions, # but the multi-strain disease model above is typed by the new `strain_ontology`. -# Crucially, we can accomplish this IF there is an appropriate morphism (map) between the type systems because post-composition by a morphism of type systems is functorial. +# Crucially, we can accomplish this IF there is an appropriate morphism (map) between the type systems +# ($P_{strain}$ and $P_{infectious}$) because post-composition by a morphism of type systems is functorial. # In this case, there is a morphism from `strain_ontology` to `infectious_ontology`, so we can form the morphism # ### Morphism from `strain_ontology` to `infectious_ontology` +# We use `oapply_typed` on a UWD representation of the `strain_ontology`, but note that we could also directly +# define the map $P_{strain}\to P_{infectious}$ with `ACSetTransformation`. + strain_ont_uwd = @relation (Uninf,Inf) where (Uninf::Pop, Inf::Pop) begin infect(Uninf, Inf, Inf, Inf) disease(Inf, Inf) @@ -226,9 +236,8 @@ to_graphviz(strain_ont_act) # To demonstrate stratification utilizing post-composition to re-type the models, we use the simple-trip geographic model. # This model is comprises a travel model and a living model. -# **Travel model between $N$ regions**\n # In this model, there are $N$ regions which people can travel between. People within the same region are able -# to infect other people in the same region. +# to infect other people in the same region. It is typed by `infectious_ontology` ($P_{infectious}$). function travel_model(n) uwd = RelationDiagram(repeat([:Pop], n)) @@ -256,8 +265,8 @@ to_graphviz(dom(travel_model(2))) # This model could itself be stratified with the SIRD model, but we want to model # persons travelling between locations while maintaining the status of where they live. -# **Living model of $N$ regions**\n -# In this model, people have the property of "Living" somewhere. +# In this model, people have the property of "Living" somewhere. +# It is also typed by `infectious_ontology` ($P_{infectious}$). function living_model(n) typed_living = pairwise_id_typed_petri(infectious_ontology, :Pop, :infect, [Symbol("Living$(i)") for i in 1:n]) @@ -266,21 +275,24 @@ end to_graphviz(dom(living_model(2))) -# **Simple trip model of $N$ regions**\n # We can stratify the living model with the travel model to get a model of someone taking a simple trip. simple_trip_model = typed_product(travel_model(2), living_model(2)) to_graphviz(dom(simple_trip_model)) -# ### Stratify SIRD-multi-strain and simple-trip models +# ### Stratify multi-strain SIRD and simple-trip models # Now, to stratify our multi-strain SIRD model with the simple-trip, we first retype the multi-strain model -# to the `infectious_ontology` by composing with the morphism we defined. +# to the `infectious_ontology` by composing with the morphism we defined. Mathematically, because +# the multi-strain SIRD model is a typed Petri net, it is a morphism $\phi : P \to P_{strain}$. +# The object `strain_ont_act` we made earlier is a morphism between $P_{strain} \to P_{infectious}$, +# so when we post-compose, we get a new morphism between $P \to P_{infectious}$. sird_strain_retyped = compose(sird_strain,strain_ont_act) -# We can now stratify. +# We can now take the typed product to get the multi-strain SIRD model stratified over the simple trip model +# of movement between different regions. sird_strain_trip = typed_product(sird_strain_retyped,simple_trip_model) @@ -288,7 +300,9 @@ to_graphviz(dom(sird_strain_trip)) # ## Define a multi-strain SIRD model with vaccination by multiple vaccine types -# We can similarly stratify the multi-strain SIRD model with the multi-vax model. +# Because the multi-vaccine model was typed according to $P_{infectious}$, our retyped +# multi-strain SIRD can be stratified according to the multiple vaccine model in a +# similar way. sird_strain_vax = typed_product(sird_strain_retyped,vax_model(2)) @@ -299,9 +313,11 @@ to_graphviz(dom(sird_strain_vax)) # If we would like to re-stratify our SIRD-strain-vax model with the simple trip model, we again face a difficulty. # Both the "vaccination" transitions of the first model and the "travel" transitions of the second # are currently typed to the `:strata` transition of the `infectious_ontology` type system. +# Naively stratifying would thus produce transitions in which persons traveled and were vaccinated simultaneously. + # To appropriately stratify, we need an additional "strata" transition to distinguish # between the two types of transitions. -# We can again use post-compostion to overcome the difficulty. +# We can again use post-compostion with the morphism between type systems to reuse our existing models. # ### Define an augmented version of the `infectious_ontology` type system with an additional "strata" transition diff --git a/docs/literate/covid/epidemiology.jl b/docs/literate/covid/epidemiology.jl index 276a960e..72e3fd27 100644 --- a/docs/literate/covid/epidemiology.jl +++ b/docs/literate/covid/epidemiology.jl @@ -17,34 +17,55 @@ using Catlab.Programs.RelationalPrograms display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75")); -# #### SIR Model: +# #### SIR Model -# define model +# We use the `@relation` macro from Catlab to create an undirected wiring diagram (UWD) +# which describes the composition syntax for the SIR model. Briefly, boxes (labeled ovals) +# represent processes which may depend on (consume or produce) resources represented by +# junctions (labeled dots). sir = @relation (s,i,r) begin infection(s,i) recovery(i,r) end display_uwd(sir) + +# To generate a Petri net (PN) from our compositional syntax, we need to apply a composition +# syntax, which assigns concrete mathematical models to each process. To compose +# with PNs, each box in the UWD will correspond to an open PN whose feet attach to +# the junctions that box is connected to. The composite PN is then constructed by +# gluing the component PNs along at the shared junctions (places). +# +# The method `oapply_epi` is used to do this composition. It is a wrapper for the `oapply` +# method, which can be used with general models, and substitutes PNs for boxes +# in the UWD according to the box name. For a list of the allowed box labels/PNs, +# please see the help file for the function. +# +# The returned object is a `MultiCospan`, where the feet are each outer port of the UWD +# and legs are morphisms which identify each outer port to a place in the composed PN. +# The apex of the multicospan is thus the composed model. #- p_sir = apex(oapply_epi(sir)) to_graphviz(p_sir) -# define initial states and transition rates, then -# create, solve, and visualize ODE problem +# Labelled vectors are used to create the initial state and reaction rate parameters for each transition. u0 = LVector(S=10, I=1, R=0); p = LVector(inf=0.4, rec=0.4); -# The C-Set representation has direct support for generating a DiffEq vector field +# The `vectorfield` method interprets the PN as describing mass-action kinetics +# with a rate constant associated to each transition, which can be used to +# simulate ODEs associated to the PN. prob = ODEProblem(vectorfield(p_sir),u0,(0.0,7.5),p); sol = solve(prob,Tsit5()) -plot(sol) +plot(sol, labels=["S" "I" "R"]) + +# #### SEIR Model -# #### SEIR Model: +# Like the SIR model, we use an undirected wiring diagram as the composition syntax, where +# junctions correspond to places and boxes to Petri nets. -# define model seir = @relation (s,e,i,r) begin exposure(s,i,e) illness(e,i) @@ -55,8 +76,8 @@ display_uwd(seir) p_seir = apex(oapply_epi(seir)) to_graphviz(p_seir) -# define initial states and transition rates, then -# create, solve, and visualize ODE problem +# Define initial states and transition rates, then +# create, solve, and visualize ODE problem: u0 = LVector(S=10, E=1, I=0, R=0); p = LVector(exp=.9, ill=.2, rec=.5); @@ -64,11 +85,12 @@ p = LVector(exp=.9, ill=.2, rec=.5); prob = ODEProblem(vectorfield(p_seir),u0,(0.0,15.0),p); sol = solve(prob,Tsit5()) -plot(sol) +plot(sol, labels=["S" "E" "I" "R"]) # #### SEIRD Model: -# define model +# We can add a deceased component and a death process to the SEIR model, specified with +# an undirected wiring diagram. seird = @relation (s,e,i,r,d) begin exposure(s,i,e) illness(e,i) @@ -80,8 +102,8 @@ display_uwd(seird) p_seird = apex(oapply_epi(seird)) to_graphviz(p_seird) -# define initial states and transition rates, then -# create, solve, and visualize ODE problem +# Define initial states and transition rates, then +# create, solve, and visualize ODE problem: u0 = LVector(S=10, E=1, I=0, R=0, D=0); p = LVector(exp=0.9, ill=0.2, rec=0.5, death=0.1); @@ -89,4 +111,4 @@ p = LVector(exp=0.9, ill=0.2, rec=0.5, death=0.1); prob = ODEProblem(vectorfield(p_seird),u0,(0.0,15.0),p); sol = solve(prob,Tsit5()) -plot(sol) +plot(sol, labels=["S" "E" "I" "R" "D"]) diff --git a/docs/literate/covid/stratification.jl b/docs/literate/covid/stratification.jl index 547d53a0..9c4975b0 100644 --- a/docs/literate/covid/stratification.jl +++ b/docs/literate/covid/stratification.jl @@ -1,4 +1,4 @@ -# # [Stratification of COVID Models](@id stratification_example) +# # [Stratification of Epidemiological Models](@id stratification_example) # #md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/covid/stratification.ipynb) @@ -8,9 +8,20 @@ using Catlab.CategoricalAlgebra using Catlab.WiringDiagrams using DisplayAs, Markdown +display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75")); + +# In this tutorial we show how to use this package to define basic epidemiological models and +# model stratification, when a submodel is replicated across a variable and allowed to interact +# in some specified way (e.g.; age-structured infection dynamics). For more details on the underlying +# theory, please see the research article ["An algebraic framework for structured epidemic modelling"](https://doi.org/10.1098/rsta.2021.0309). + # ## Define basic epidemiology model -# Define the type system for infectious disease models +# One of the fundamental parts of getting model stratification to work is defining a "type system" +# for the models one is interested in, which ensures that spurious or erroneous transitions are +# not generated when we construct the stratified model (e.g.; no vector to vector transmission in a malaria model). +# In the research article this tutorial follows, the Petri net that defines the type system is ``P_{type}``. +# We define ``P_{infectious}`` below. const infectious_ontology = LabelledPetriNet( [:Pop], @@ -21,9 +32,23 @@ const infectious_ontology = LabelledPetriNet( to_graphviz(infectious_ontology) -# Define a simple SIRD model with reflexive transitions typed as `:strata` to indicate which states can be stratified +# More details are given in the article, but briefly, ``P_{infectious}`` has a single place, which represents +# a potentially stratified population, and 3 "types" of transitions: "infect" is for infections, "disease" for changes +# in disease state (e.g.; from E to I), and "strata" for changes not in disease state (e.g.; from region ``i`` to region ``j``). +# +# A typed Petri net is a morphism ``\phi: P \to P_{type}``. Here we define a Petri net representing the SIRD model +# typed over ``P_{infectious}``. We start by defining an undirected wiring diagram (UWD) which represents the processes +# and states in the SIRD model; note that the name of each box corresponds to one of the transitions in ``P_{infectious}``, +# and that we specify the types of the junctions as `Pop` in the `where` clause. +# +# Next we use `oapply_typed` which produces a typed Petri net by composing transitions based on ``P_{infectious}``. +# The typed Petri net is modified with reflexive transitions typed as `:strata` to indicate which states can be stratified # Here we add reflexive transitions to the susceptible, infected, and recovered populations but we leave out the dead -# population because they cannote do things such as get vaccinated or travel between regions. +# population as no individuals entering that compartment may leave. +# +# Calling `to_graphviz` on the resulting object will display $\phi : P \to P_{infectious}$. We could apply `dom` (domain) +# or `codom` (codomain) on the object to extract $P$ or $P_{infectious}$, respectively, if we just wanted to plot +# a single Petri net. sird_uwd = @relation (S,I,R,D) where (S::Pop, I::Pop, R::Pop, D::Pop) begin infect(S, I, I, I) @@ -34,29 +59,46 @@ end sird_model = oapply_typed(infectious_ontology, sird_uwd, [:infection, :recovery, :death]) sird_model = add_reflexives(sird_model, [[:strata], [:strata], [:strata], []], infectious_ontology) -to_graphviz(dom(sird_model)) +to_graphviz(sird_model) # ## Define intervention models # ### Masking model +# Let's say we want to stratify the populations in the SIRD model by masking or non-masking; or, generally +# by any behavioral change which is reversible and does not preclude the possibility of infection in either +# state. To generate our stratification scheme, we begin with a UWD describing the strata levels we are +# interested in (Masked and UnMasked), and the transitions possible in each level. Note that we give +# transitions of type `strata` which allow masking or unmasking. Unmasked individuals may transmit +# disease to either masked or unmasked individuals, but masked persons cannot transmit disease. +# +# We again use `oapply_typed` to generate the stratification scheme ``\phi^{'}:P^{'}\to P_{infectious}``, +# adding reflexive `disease` transitions to the typed Petri net as changes in disease status may occur +# in either stratum. + masking_uwd = @relation (M,UM) where (M::Pop, UM::Pop) begin - disease(M, UM) - disease(UM, M) + strata(M, UM) + strata(UM, M) infect(M, UM, M, UM) infect(UM, UM, UM, UM) end mask_model = oapply_typed(infectious_ontology, masking_uwd, [:unmask, :mask, :infect_um, :infect_uu]) -mask_model = add_reflexives(mask_model, [[:strata], [:strata]], infectious_ontology) +mask_model = add_reflexives(mask_model, [[:disease], [:disease]], infectious_ontology) to_graphviz(dom(mask_model)) +# To generate the stratified SIRD model, we use `typed_product`. As described in the article, a stratification of +# a model can be seen as a pullback of ``\phi`` and ``\phi^{'}``, or a product in the slice category ``Petri/P_{type}``. # Stratify our SIRD model on this masking model to get a model of SIRD with masking: typed_product(sird_model, mask_model) |> dom |> to_graphviz # ### Vaccine model +# By changing ``P^{'}`` slightly, we can generate a stratification scheme to model vaccination. In this case, +# the transition between strata is irreversible (i.e.; the vaccine is not leaky), infection may occur between any +# pair of individuals, and change in disease state may occur in any strata. + vax_uwd = @relation (UV,V) where (UV::Pop, V::Pop) begin strata(UV, V) infect(V, V, V, V) @@ -69,7 +111,7 @@ vax_model = add_reflexives(vax_model, [[:disease], [:disease]], infectious_ontol to_graphviz(dom(vax_model)) -# Stratify our SIRD model on this vaccine model to get a model of SIRD with a vaccination rate: +# Once again, the typed product, or pullback gives the SIRD model stratified by vaccination. typed_product(sird_model, vax_model) |> dom |> to_graphviz @@ -108,9 +150,10 @@ typed_product(sird_model, mask_vax_model) |> dom |> to_graphviz # ### Travel model between $N$ regions -# For this model we can use a julia function to programmatically build up our undirected wiring diagram for defining this model. -# Here we want there to be $N$ regions in which people can travel between each region and people within the same region are able -# to infect other people in the same region. +# In many cases, it is not practical to construct by hand the undirected wiring diagram used for composition. Here we demonstrate +# how to use the imperative interface provided by Catlab to construct a UWD describing a model of a population spread out amongst $N$ +# regions. Individuals can travel between regions, but disease transmission is only possible between individuals in the same region. +# The result is a typed Petri net describing this travel model. function travel_model(n) uwd = RelationDiagram(repeat([:Pop], n)) @@ -144,6 +187,9 @@ typed_product(sird_model, travel_model(2)) |> dom |> to_graphviz # For this model we can use a julia function to programmatically build up our model where people have the property of living somewhere # and we are modelling them travelling between locations while maintaining the status of where they live. Here we can actually just # define the model of having a "Living" status and stratify it with the previously defined travel model to get a model of someone taking a simple trip. +# In the "living model" we allow for infections between individuals who live at different places; when we take the typed +# product with the simple trip model this will result in a model where infection is allowed between individuals at +# the same region, regardless of their home region. function living_model(n) typed_living = pairwise_id_typed_petri(infectious_ontology, :Pop, :infect, [Symbol("Living$(i)") for i in 1:n]) diff --git a/docs/literate/predation/lotka-volterra.jl b/docs/literate/predation/lotka-volterra.jl index 7f1ec6a1..aeb389ee 100644 --- a/docs/literate/predation/lotka-volterra.jl +++ b/docs/literate/predation/lotka-volterra.jl @@ -15,8 +15,21 @@ using Catlab.Programs.RelationalPrograms display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75")); -# #### Step 1: Define the building block Petri nets needed to construct the model +# #### Introduction +# In this tutorial we use the [Lotka-Volterra model](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations), a classic +# mathematical model from ecology, to demonstrate basic conceps of compositional model building using AlgebraicPetri. + +# #### Step 1: Define the building block Petri nets needed to construct the model. +# +# These are the basic Petri nets which will be substituted as concrete models in the compositional syntax. +# There are three processes which we need to represent explicitly: birth, predation, and death. +# Because these models will be composed, we use `Open` to generate an open Petri net. When given only a single argument, +# `Open` generates a structured multicospan, which is an object containing the original Petri net ``S`` (accessed by `apex`), +# a list of finite sets ``A_{1},\dots,A_{n}`` (accessed by `feet`), and morphisms ``A_{1}\to S,\dots,A_{n}\to S`` (accessed by `legs`). +# The feet are where the Petri net may interact with other systems, justifying the term "open" Petri nets. +# +# The basic Petri nets are generated and the apex of each open Petri net is displayed. birth_petri = Open(PetriNet(1, 1=>(1,1))); to_graphviz(birth_petri) #- @@ -26,8 +39,18 @@ to_graphviz(predation_petri) death_petri = Open(PetriNet(1, 1=>())); to_graphviz(death_petri) +# #### Step 2: Specify interaction between submodels using a relational syntax -# #### Step 2: Generate models using a relational syntax +# To stitch independent open Petri nets together into a combined model, we need to specify an interaction pattern (composition syntax). +# This is specified as an undirected wiring diagram (UWD) using the `@relation` macro from Catlab.jl, please see the documentation +# there for more information. +# +# The UWD produced here consists of 3 "boxes" and 2 "junctions". Boxes connect to junctions via wires based on the function-like syntax +# in the macro. Wires that lead off the display are "outer ports" indicating how output may be collected from the combined system. +# +# We use `oapply` to substitute the open Petri nets constructed earlier into each box of the UWD, where the feet of each open Petri net +# are identified with the wires from each box to shared junctions. The result of `oapply` is, once again, a multicospan, where the legs +# now give the map between outer ports of the composed model to places. lotka_volterra = @relation (wolves, rabbits) begin birth(rabbits) @@ -40,16 +63,19 @@ lv_dict = Dict(:birth=>birth_petri, :predation=>predation_petri, :death=>death_p lotka_petri = apex(oapply(lotka_volterra, lv_dict)) to_graphviz(lotka_petri) -# Generate appropriate vector fields, define parameters, and visualize solution +# We can now define an initial state `u0`, reaction rate constants for each transition `p`, and use the method `vectorfield` +# to create a function we can pass to a differential equation solver to simulate and plot a trajectory of our +# composed Lotka-Volterra model. u0 = [100, 10]; p = [.3, .015, .7]; prob = ODEProblem(vectorfield(lotka_petri),u0,(0.0,100.0),p); sol = solve(prob,Tsit5(),abstol=1e-8); -plot(sol) +plot(sol, labels=["Rabbits" "Wolves"]) # #### Step 3: Extend your model to handle more complex phenomena -# such as a small food chain between little fish, big fish, and sharks + +# By making a larger UWD to describe, we can describe a small food chain between little fish, big fish, and sharks. dual_lv = @relation (fish, Fish, Shark) begin birth(fish) @@ -69,4 +95,59 @@ u0 = [100, 10, 2]; p = [.3, .015, .7, .017, .35]; prob = ODEProblem(vectorfield(dual_lv_petri),u0,(0.0,100.0),p); sol = solve(prob,Tsit5(),abstol=1e-6); -plot(sol) +plot(sol, label=["Little fish" "Big fish" "Sharks"]) + +# However, at this point it is natural to be concerned about the scalability of needing to define +# the UWD for larger models; is there a way to compose UWDs themselves to get a larger UWD? +# As you may have guessed, the answer is yes. Let's see how. +# +# Let's say we want to include another predator of rabbits into the model, hawks. Hawks +# will behave like flying wolves in this setup, but they don't need to. + +lotka_volterra_hawk = @relation (hawks, rabbits) begin + birth(rabbits) + predation(rabbits, hawks) + death(hawks) +end +display_uwd(lotka_volterra_hawk) + +# We would like to "glue" the two UWDs describing terrestrial and aerial predation together +# along the rabbits junction and birth box. We cannot simply take a disjoint union of the +# two UWDs (coproduct) as then we would be modelling two disconnected systems. Instead we +# need to perform a categorical operation called a "pushout", which is like a disjoint union +# but where certain elements are glued together. +# +# We specify this overlap between the two systems where we will glue them together below. + +rabbit_uwd = @relation (rabbits,) begin + birth(rabbits) +end +display_uwd(rabbit_uwd) + +# Next we need to define mappings from this overlap into each of the UWDs we want to glue together. +# We use the `ACSetTransformation` method from Catlab for this, which defines a natural transformation +# between two C-Sets. The result is a "span", where the apex is the overlap, and the acset transformations +# are arrows from the apex into each of the original UWDs. We take the pushout of the span +# to glue the UWDs together along the apex. + +hawk_transform = ACSetTransformation((Box=[1], Junction=[2], Port=[1], OuterPort=[2]), rabbit_uwd, lotka_volterra_hawk) +wolf_transform = ACSetTransformation((Box=[1], Junction=[2], Port=[1], OuterPort=[2]), rabbit_uwd, lotka_volterra) +lotka_volterra_composed = ob(pushout(hawk_transform, wolf_transform)) + +display_uwd(lotka_volterra_composed) + +# Now that we have the composition syntax for the combined terrestrial-aerial predation system, we +# can substitute concrete mathematical models (open Petri nets) into each box as before to +# produce a concrete model of the composed system. + +lotka_petri = apex(oapply(lotka_volterra_composed, lv_dict)) +to_graphviz(lotka_petri) + +# We may now again interpret the model as a differential equation, +# solve, and plot the solution. + +u0 = [10.0, 100.0, 50.0]; +p = [.3, .015, .7, 0.01, 0.5]; +prob = ODEProblem(vectorfield(lotka_petri),u0,(0.0,100.0),p); +sol = solve(prob,Tsit5()); +plot(sol, label=["Rabbits" "Wolves" "Hawks"]) diff --git a/src/Epidemiology.jl b/src/Epidemiology.jl index e1224753..2554bc70 100644 --- a/src/Epidemiology.jl +++ b/src/Epidemiology.jl @@ -49,7 +49,7 @@ epi_dict = Dict(:infection=>infection, :exposure=>exposure, :illness=>illness, : """ oapply_epi(ex, args...) Generates a LabelledPetriNet under a composition pattern described by the -undirected wiring diagram `ex`. This requires that the nodes in `ex` are only +undirected wiring diagram `ex`. This requires that the boxes in `ex` are only labelled with labels from the following set: ``` [:infection, :exposure, :illness, :recovery, :death] diff --git a/src/TypedPetri.jl b/src/TypedPetri.jl index 83ba86c1..6791de61 100644 --- a/src/TypedPetri.jl +++ b/src/TypedPetri.jl @@ -61,6 +61,8 @@ where each of the boxes is labeled by a symbol that matches the label of a transition in the petri net. Then produces a petri net given by colimiting the transitions together, and returns the ACSetTransformation from that Petri net to the type system. + +The `tnames` argument renames the transitions in the resulting typed Petri net. """ function oapply_typed(type_system::LabelledPetriNet, uwd, tnames::Vector{Symbol}) if junction(uwd, outer=true) != junctions(uwd) @@ -97,7 +99,7 @@ function oapply_typed(type_system::LabelledPetriNet, uwd, tnames::Vector{Symbol} ) end -r""" +""" Modify a typed petri net to add "reflexive transitions". These are transitions which go from one species to itself, so they don't change the mass action semantics, but they are important for stratification. The idea behind this is similar to the fact that the product of the graph *----* with itself is