diff --git a/episodes/superspreading-simulate.Rmd b/episodes/superspreading-simulate.Rmd
index 8f2ce5e8..a99f3762 100644
--- a/episodes/superspreading-simulate.Rmd
+++ b/episodes/superspreading-simulate.Rmd
@@ -705,14 +705,27 @@ aggregate_chains <- multiple_epichains %>%
   dplyr::as_tibble() %>%
   # get the round number (day) of infection times
   dplyr::mutate(day = ceiling(time)) %>%
-  # count the daily number of cases in each chain
-  dplyr::count(chain, day, name = "cases") %>%
+  # count the number of daily incident cases in each chain
+  dplyr::count(chain, day, name = "cases_new") %>%
   # calculate the cumulative number of cases for each chain
   dplyr::group_by(chain) %>%
-  dplyr::mutate(cases_cumsum = cumsum(cases)) %>%
+  dplyr::mutate(cases_cumsum = cumsum(cases_new)) %>%
   dplyr::ungroup()
 ```
 
+To better understand the output we just produced, let's print more lines of it. The relevant columns are:
+
+- **`days`**: the discretized time step of the simulation (continuous time binned into days)
+- **`cases_new`**: the number of new incident cases on that day
+- **`cases_cumsum`**: the cumulative sum of cases per day, computed within each transmission chain
+
+```r
+# Print a couple of chains to better understand the output
+aggregate_chains %>% 
+  arrange(chain) %>% 
+  print(n=50)
+```
+
 Before the plot, let's create a summary table with the total time duration and size of each chain. We can use the `{dplyr}` "combo" of `group_by()`, `summarise()` and `ungroup()`:
 
 ```{r}
@@ -788,7 +801,6 @@ This output should give us equivalent results to `summary(multiple_epichains)`:
 ```{r}
 # number of chains that reached the 100-case threshold
 summary_chains %>%
-  dplyr::arrange(desc(cases_total)) %>%
   dplyr::filter(cases_total > 100)
 ```