Skip to content

Commit

Permalink
Merge pull request #169 from slwu89/docs
Browse files Browse the repository at this point in the history
Improvements to documentation
  • Loading branch information
jpfairbanks authored Mar 20, 2024
2 parents cd46566 + 76fa519 commit 8828177
Show file tree
Hide file tree
Showing 6 changed files with 226 additions and 59 deletions.
64 changes: 40 additions & 24 deletions docs/literate/covid/disease_strains.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand All @@ -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)
Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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))
Expand Down Expand Up @@ -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])
Expand All @@ -266,29 +275,34 @@ 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)

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))

Expand All @@ -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

Expand Down
52 changes: 37 additions & 15 deletions docs/literate/covid/epidemiology.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -55,20 +76,21 @@ 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);

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)
Expand All @@ -80,13 +102,13 @@ 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);

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"])
Loading

0 comments on commit 8828177

Please sign in to comment.