diff --git a/inst/starter_models/sir_age/README.Rmd b/inst/starter_models/sir_age/README.Rmd
new file mode 100644
index 000000000..209cba5a2
--- /dev/null
+++ b/inst/starter_models/sir_age/README.Rmd
@@ -0,0 +1,168 @@
+---
+title: "Age-stratified SIR"
+index_entry: "An age-stratified SIR model"
+bibliography: ../../references.bib
+link-citations: TRUE
+author: Irena Papst, Steve Walker
+output: 
+  github_document:
+    toc: true
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  warning = FALSE,
+  message = FALSE,
+  comment = "#>",
+  fig.path = "./figures/"
+)
+```
+
+This is an expansion of the [sir](https://github.com/canmod/macpan2/tree/main/inst/starter_models/sir) model to include age-stratification; each state variable is subdivided by age, and transmission is governed by a age-based contact matrix, as described in [@Mistry2021].
+
+# Packages Used and Settings
+
+The code in this article uses the following packages.
+
+```{r packages, message=FALSE, warning=FALSE}
+library(ggplot2)
+library(dplyr)
+library(forcats)
+library(scales)
+library(macpan2)
+```
+
+# Model Specification
+
+This model has been specified in the `sir_age` directory [here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_age/tmb.R) and is accessible from the `macpan2` model library (see [Example Models](https://canmod.github.io/macpan2/articles/example_models.html) for details). We can read in the model specification using the `mp_tmb_library` command.
+```{r model_spec}
+spec = mp_tmb_library(
+    "starter_models"
+  , "sir_age"
+  , package = "macpan2"
+)
+```
+
+This specification can be used to draw the following flow diagram using code found in the [source for this article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_age/README.Rmd).
+
+```{r diagram, echo = FALSE, eval = FALSE, fig.height = 2, fig.width = 8}
+system.file("utilss", "box-drawing.R", package = "macpan2") |> source()
+layout = mp_layout_paths(spec)
+(layout
+  |> plot_flow_diagram(show_flow_names = TRUE)
+  |> draw_inflows(layout, show_labels = TRUE)
+  |> draw_outflows(layout, show_labels = TRUE, pattern_mutate = "^([a-z]*)(_.*)", pattern_replace = "\\1")
+)
+```
+
+
+# States
+
+| variable | description                       |
+| -------- | --------------------------------- |
+| S        | Number of susceptible individuals |
+| I        | Number of infectious individuals  |
+| R        | Number of recovered individuals   |
+
+Each state variable is a vector with one entry per age group. 
+
+The size of the population in each age groups is given by the components of $N = S + I + R$.
+
+# Parameters
+
+| variable | description                  |
+| -------- | ---------------------------- |
+| $\tau$  | transmissibility, _i.e._, the proportion of contacts with infecteds that yield infection |
+| $\gamma$ | per capita recovery rate     |
+| $M$    | matrix of contact rates between age groups, _e.g._ $M_{ij}$ is the rate of contact between susceptibles of age $i$ and infecteds of age $j$        |
+
+$\tau$ and $\gamma$ can be specified by age (as vectors), though in the example below we use a scalar value that applies to all age groups.
+
+# Dynamics 
+
+$$
+\begin{align*}
+\frac{dS}{dt} &= - \lambda S \\
+\frac{dI}{dt} &= \lambda S - \gamma I \\
+\frac{dR}{dt} &= \gamma I
+\end{align*}
+$$
+
+Here $\lambda = \tau M \frac{I}{N}$.
+
+# Simulation Example
+
+```{r plot_sim}
+# set numerical values
+N <- rep(1e5, 3) # total population by age
+I0 <- 100 # initial infecteds
+
+M <- matrix(c(
+  10, 1, 1, 
+  1, 5, 1, 
+  1, 1, 2
+  ), nrow = 3) # contact matrix in units of avg number of contacts/person/day (symmetric contact rates because we will pick two age groups of the same size)
+gamma <- 1/5 # 5 days avg infectious period; same gamma across age groups
+
+R0 <- 2
+# spectral radius of contact matrix for transmissibility calculation
+rho <- max(abs(eigen(M)$values))
+# calculate transmissibility given R0, gamma, M
+tau <- R0*gamma/rho # transmissibility (proportion of contacts with infectious individuals that yield infection)
+
+# variable labellers
+matrix_labeller <-  c(
+  "S" = "Susceptible",
+  "I" = "Infectious",
+  "R" = "Recovered"
+)
+row_labeller <- c(
+  "0" = "Children (ages 0-17)",
+  "1" = "Adults (ages 18-64)",
+  "2" = "Seniors (ages 65+)"
+)
+
+# simulate & visualize
+(spec
+  # simulate
+  |> mp_simulator(
+    time_steps = 150
+  , outputs = mp_state_vars(spec)
+  )
+  |> mp_trajectory()
+  # enforce natural order of variables for plotting
+  |> dplyr::mutate(
+    matrix = factor(matrix, levels = mp_state_vars(spec))
+   , row = forcats::as_factor(row)
+   )
+  # visualize trajectory
+  |> ggplot2::ggplot(
+    mapping = ggplot2::aes(x = time, y = value, colour = row)
+  ) 
+  + ggplot2::geom_line() +
+  + ggplot2::facet_wrap(
+      ~ matrix, 
+      labeller = ggplot2::labeller(
+        matrix = matrix_labeller
+      ),
+      scales = "free_y")
+  + ggplot2::scale_colour_discrete(
+      labels = row_labeller
+    )
+  + ggplot2::labs(
+      title = "Simulation of the age-structured SIR model with two age groups",
+      subtitle = paste("Infection is seeded with", I0, "infectious children in a total population of", scales::label_number(scale_cut = scales::cut_short_scale())(sum(N))),
+      x = "Days from beginning of outbreak",
+      y = "Number of individuals"
+    )
+  + ggplot2::theme(
+      legend.position = "bottom"
+    , legend.title = ggplot2::element_blank()
+    )
+)
+```
+
+# References
diff --git a/inst/starter_models/sir_age/references.bib b/inst/starter_models/sir_age/references.bib
new file mode 100644
index 000000000..ce5f6f045
--- /dev/null
+++ b/inst/starter_models/sir_age/references.bib
@@ -0,0 +1,13 @@
+@article{Mistry2021,
+  title = {Inferring high-resolution human mixing patterns for disease modeling},
+  volume = {12},
+  ISSN = {2041-1723},
+  url = {http://dx.doi.org/10.1038/s41467-020-20544-y},
+  DOI = {10.1038/s41467-020-20544-y},
+  number = {1},
+  journal = {Nature Communications},
+  publisher = {Springer Science and Business Media LLC},
+  author = {Mistry,  Dina and Litvinova,  Maria and Pastore y Piontti,  Ana and Chinazzi,  Matteo and Fumanelli,  Laura and Gomes,  Marcelo F. C. and Haque,  Syed A. and Liu,  Quan-Hui and Mu,  Kunpeng and Xiong,  Xinyue and Halloran,  M. Elizabeth and Longini,  Ira M. and Merler,  Stefano and Ajelli,  Marco and Vespignani,  Alessandro},
+  year = {2021},
+  month = jan 
+}
\ No newline at end of file
diff --git a/inst/starter_models/sir_age/tmb.R b/inst/starter_models/sir_age/tmb.R
new file mode 100644
index 000000000..280186fc1
--- /dev/null
+++ b/inst/starter_models/sir_age/tmb.R
@@ -0,0 +1,35 @@
+library(macpan2)
+
+# SIR with age, where age structure only affects transmission, as 
+# defined here: https://www.nature.com/articles/s41467-020-20544-y#Sec7
+
+# define the model
+spec = mp_tmb_model_spec(
+  before = list(
+    S ~ N - I - R # calculate population size before the simulation loop begins to avoid having to specify a value for it by hand in the defaults list
+  ),
+  during = list(
+    # force of infection
+    lambda ~ tau * M %*% (I / N),
+    # infection
+    mp_per_capita_flow(
+      from = "S", to = "I",
+      rate = "lambda",
+      abs_rate = "infection"
+    ),
+    # recovery
+    mp_per_capita_flow(
+      from = "I", to = "R",
+      rate = "gamma",
+      abs_rate = "recovery"
+    )
+  ),
+  default = list(
+    tau = tau,
+    M = M,
+    gamma = gamma, # 5 days avg infectious period; same gamma across age groups 
+    N = N,
+    I = c(I0, 0, 0), # seed infection in age group 1
+    R = rep(0, 3)
+  )
+)