Introduction

One question that has come up a lot this semester is, “When do you need to include a latent (\(E\)) compartment in your dynamical model?”

The short answer is when the duration of the latent period about the same magnitude as the duration of the infectious period or when it is a “significant” amount of time. What that actually means really depends on your question and your model. This notebook is just a brief but systematic look at how changing your latent period affects the dynamics of a disease.

The basic SEIR model

For this, let’s assume there are no deaths, no births, and lifelong immunity. That is, our model can be written as:

\[ \begin{aligned} \frac{dS}{dt} &= - \beta S I \\ \frac{dE}{dt} &= \beta S I - \alpha E \\ \frac{dI}{dt} &= \alpha E - \gamma I \\ \frac{dR}{dt} &= \gamma I \end{aligned} \]

Like all models in this class, assume the parameters are exponentially distributed so the average time in any compartment is the inverse of the rate out of the compartment. That is, the duration of incubation is \(\frac{1}{\alpha}\) and duration of infectiousness is \(\frac{1}{\gamma}\).

Simulation set up

Given the above, let’s see what happens we vary \(\alpha\) and \(\gamma\) while fixing \(\beta\) and expressing all compartments as proportions.

Define our model in R

library(deSolve)
library(tidyverse)
SEIR <- function(t, x, parms){
    with(as.list(c(parms, x)), {
        
        dS <- -(beta * S * I) 
        dE <- +(beta * S * I) - (a * E)
        dI <- +(a * E) - (g * I)
        dR <- +(g * I)
        
        der <- c(dS, dE, dI, dR)
        
        return(list(der))
    }) 
}

Duration of incubation and infectious periods

Let’s try different length infectious periods. Specifically, let’s set \(\alpha \in \{0.05, 0.1, 0.2, 0.5, 0.99\}\), which corresponds to durations of infectiousness of \(\{20, 10, 5, 2, 1\}\) time units. For each infectious period, we will also vary the duration of the incubation period as a ratio of the infectious period. For example, let’s set the ratios of the incubation period to be \(\{0.0001, 0.1, 0.5, 1, 2\}\), which means for an infectious period of 20 days (\(\alpha = 0.05\)), we will run the simulation with incubation periods of \(\{\approx0, 2, 10, 20, 40\}\) days.

alphas <- c(.05, .1, .2, .5, .99)
gamma_props = c(.0001, .1, .5, 1, 2)
dt <- seq(0, 1000, 1)
inits <- c(S = (10^6 - 1)/10^6, 
           E = 0, 
           I = 1/10^6, 
           R = 0)
holder <- NULL
for (alpha in alphas) {
    for (g_prop in gamma_props) {
        gamma <- g_prop * alpha
        
        parms <- c(beta = 3, a = alpha, g = gamma)
        
        simulation <- as.data.frame(ode(y = inits, times = dt, 
                                func = SEIR, parms = parms))
        
        holder <- rbind(holder, cbind(t = simulation[, "time"], 
                                      infected = simulation[, "I"], 
                                      alpha = alpha, g_prop = g_prop))
    }
    
}

We now have an object which contains the infectious curve for all 25 combinations of infectious period and incubation periods.

We can plot it like so:

holder <- holder %>% 
    as_tibble() %>% 
    mutate(infectious_duration = round(1/alpha, 0))
ggplot(holder, aes(x = t, y = infected)) + 
    geom_line() + 
    facet_grid(infectious_duration ~ g_prop) + 
    scale_x_continuous(limits = c(0, 100)) + 
    theme_minimal() + 
    labs(x = "Time", y = "Proportion Infected", 
         title = "Infectious curves for various infectious and incubation periods")

Each panel represents some instance of \(\alpha\) and \(\gamma\) from our sets. Each row is one set of infectious periods (i.e., the top is an infecitous period of 1 day, the bottom is 20 days). The columns show the changing incubation periods as a ratios of the infectious periods.

If you focus on any row and go from left to right, you can see that as the latent period gets longer (relative to the infectious period), the epidemic peak (maximum number of infected) gets lower and a little more delayed.

If you focus on any column and go top to bottom, you can see that as you increase the duration of infectiousness, the epidemic peak is shifted to the right (because we are also increasing the amount of time they stay in the latent period), and the “peak” is generally fatter because the time spent in the compartment is longer. It does not, however, seem like there is any impact on the height of the peak.

Take home

Latent periods matter as they become longer relative to your infectious period. For latent periods lasting \(\lt 10\%\) of the infectious period, modelers will often leave it out for simplicity. You should almost certainly include it if it is greater than that though.

LS0tCnRpdGxlOiAiVW5kZXJzdGFuZGluZyB0aGUgaW1wYWN0IG9mIGxhdGVudCBwZXJpb2RzIGluIFNFSVIgbW9kZWxzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIyBJbnRyb2R1Y3Rpb24KCk9uZSBxdWVzdGlvbiB0aGF0IGhhcyBjb21lIHVwIGEgbG90IHRoaXMgc2VtZXN0ZXIgaXMsICIqV2hlbiogZG8geW91IG5lZWQgdG8gaW5jbHVkZSBhIGxhdGVudCAoJEUkKSBjb21wYXJ0bWVudCBpbiB5b3VyIGR5bmFtaWNhbCBtb2RlbD8iCgpUaGUgc2hvcnQgYW5zd2VyIGlzIHdoZW4gdGhlIGR1cmF0aW9uIG9mIHRoZSBsYXRlbnQgcGVyaW9kIGFib3V0IHRoZSBzYW1lIG1hZ25pdHVkZSBhcyB0aGUgZHVyYXRpb24gb2YgdGhlIGluZmVjdGlvdXMgcGVyaW9kIG9yIHdoZW4gaXQgaXMgYSAic2lnbmlmaWNhbnQiIGFtb3VudCBvZiB0aW1lLiBXaGF0IHRoYXQgYWN0dWFsbHkgbWVhbnMgcmVhbGx5IGRlcGVuZHMgb24geW91ciBxdWVzdGlvbiBhbmQgeW91ciBtb2RlbC4gVGhpcyBub3RlYm9vayBpcyBqdXN0IGEgYnJpZWYgYnV0IHN5c3RlbWF0aWMgbG9vayBhdCBob3cgY2hhbmdpbmcgeW91ciBsYXRlbnQgcGVyaW9kIGFmZmVjdHMgdGhlIGR5bmFtaWNzIG9mIGEgZGlzZWFzZS4KCiMjIFRoZSBiYXNpYyBTRUlSIG1vZGVsCgpGb3IgdGhpcywgbGV0J3MgYXNzdW1lIHRoZXJlIGFyZSBubyBkZWF0aHMsIG5vIGJpcnRocywgYW5kIGxpZmVsb25nIGltbXVuaXR5LiBUaGF0IGlzLCBvdXIgbW9kZWwgY2FuIGJlIHdyaXR0ZW4gYXM6CgokJApcYmVnaW57YWxpZ25lZH0KXGZyYWN7ZFN9e2R0fSAmPSAtIFxiZXRhIFMgSSBcXApcZnJhY3tkRX17ZHR9ICY9IFxiZXRhIFMgSSAtIFxhbHBoYSBFIFxcClxmcmFje2RJfXtkdH0gJj0gXGFscGhhIEUgLSBcZ2FtbWEgSSBcXApcZnJhY3tkUn17ZHR9ICY9IFxnYW1tYSBJClxlbmR7YWxpZ25lZH0KJCQKCkxpa2UgYWxsIG1vZGVscyBpbiB0aGlzIGNsYXNzLCBhc3N1bWUgdGhlIHBhcmFtZXRlcnMgYXJlIGV4cG9uZW50aWFsbHkgZGlzdHJpYnV0ZWQgc28gdGhlIGF2ZXJhZ2UgdGltZSBpbiBhbnkgY29tcGFydG1lbnQgaXMgdGhlIGludmVyc2Ugb2YgdGhlIHJhdGUgb3V0IG9mIHRoZSBjb21wYXJ0bWVudC4gVGhhdCBpcywgdGhlIGR1cmF0aW9uIG9mIGluY3ViYXRpb24gaXMgJFxmcmFjezF9e1xhbHBoYX0kIGFuZCBkdXJhdGlvbiBvZiBpbmZlY3Rpb3VzbmVzcyBpcyAkXGZyYWN7MX17XGdhbW1hfSQuCgojIyBTaW11bGF0aW9uIHNldCB1cAoKR2l2ZW4gdGhlIGFib3ZlLCBsZXQncyBzZWUgd2hhdCBoYXBwZW5zIHdlIHZhcnkgJFxhbHBoYSQgYW5kICRcZ2FtbWEkIHdoaWxlIGZpeGluZyAkXGJldGEkIGFuZCBleHByZXNzaW5nIGFsbCBjb21wYXJ0bWVudHMgYXMgcHJvcG9ydGlvbnMuCgojIyMgRGVmaW5lIG91ciBtb2RlbCBpbiBgUmAKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCByZXN1bHRzPSJoaWRlIn0KbGlicmFyeShkZVNvbHZlKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKClNFSVIgPC0gZnVuY3Rpb24odCwgeCwgcGFybXMpewogICAgd2l0aChhcy5saXN0KGMocGFybXMsIHgpKSwgewogICAgICAgIAogICAgICAgIGRTIDwtIC0oYmV0YSAqIFMgKiBJKSAKICAgICAgICBkRSA8LSArKGJldGEgKiBTICogSSkgLSAoYSAqIEUpCiAgICAgICAgZEkgPC0gKyhhICogRSkgLSAoZyAqIEkpCiAgICAgICAgZFIgPC0gKyhnICogSSkKICAgICAgICAKICAgICAgICBkZXIgPC0gYyhkUywgZEUsIGRJLCBkUikKICAgICAgICAKICAgICAgICByZXR1cm4obGlzdChkZXIpKQogICAgfSkgCn0KYGBgCgoKIyMjIER1cmF0aW9uIG9mIGluY3ViYXRpb24gYW5kIGluZmVjdGlvdXMgcGVyaW9kcwoKTGV0J3MgdHJ5IGRpZmZlcmVudCBsZW5ndGggaW5mZWN0aW91cyBwZXJpb2RzLiBTcGVjaWZpY2FsbHksIGxldCdzIHNldCAkXGFscGhhIFxpbiBcezAuMDUsIDAuMSwgMC4yLCAwLjUsIDAuOTlcfSQsIHdoaWNoIGNvcnJlc3BvbmRzIHRvIGR1cmF0aW9ucyBvZiBpbmZlY3Rpb3VzbmVzcyBvZiAkXHsyMCwgMTAsIDUsIDIsIDFcfSQgdGltZSB1bml0cy4gRm9yIGVhY2ggaW5mZWN0aW91cyBwZXJpb2QsIHdlIHdpbGwgYWxzbyB2YXJ5IHRoZSBkdXJhdGlvbiBvZiB0aGUgaW5jdWJhdGlvbiBwZXJpb2QgYXMgYSByYXRpbyBvZiB0aGUgaW5mZWN0aW91cyBwZXJpb2QuIEZvciBleGFtcGxlLCBsZXQncyBzZXQgdGhlIHJhdGlvcyBvZiB0aGUgaW5jdWJhdGlvbiBwZXJpb2QgdG8gYmUgJFx7MC4wMDAxLCAwLjEsIDAuNSwgMSwgMlx9JCwgd2hpY2ggbWVhbnMgZm9yIGFuIGluZmVjdGlvdXMgcGVyaW9kIG9mIDIwIGRheXMgKCRcYWxwaGEgPSAwLjA1JCksIHdlIHdpbGwgcnVuIHRoZSBzaW11bGF0aW9uIHdpdGggaW5jdWJhdGlvbiBwZXJpb2RzIG9mICRce1xhcHByb3gwLCAyLCAxMCwgMjAsIDQwXH0kIGRheXMuCgpgYGB7cn0KYWxwaGFzIDwtIGMoLjA1LCAuMSwgLjIsIC41LCAuOTkpCmdhbW1hX3Byb3BzID0gYyguMDAwMSwgLjEsIC41LCAxLCAyKQpgYGAKCgpgYGB7cn0KZHQgPC0gc2VxKDAsIDEwMDAsIDEpCmluaXRzIDwtIGMoUyA9ICgxMF42IC0gMSkvMTBeNiwgCiAgICAgICAgICAgRSA9IDAsIAogICAgICAgICAgIEkgPSAxLzEwXjYsIAogICAgICAgICAgIFIgPSAwKQoKaG9sZGVyIDwtIE5VTEwKCmZvciAoYWxwaGEgaW4gYWxwaGFzKSB7CiAgICBmb3IgKGdfcHJvcCBpbiBnYW1tYV9wcm9wcykgewogICAgICAgIGdhbW1hIDwtIGdfcHJvcCAqIGFscGhhCiAgICAgICAgCiAgICAgICAgcGFybXMgPC0gYyhiZXRhID0gMywgYSA9IGFscGhhLCBnID0gZ2FtbWEpCiAgICAgICAgCiAgICAgICAgc2ltdWxhdGlvbiA8LSBhcy5kYXRhLmZyYW1lKG9kZSh5ID0gaW5pdHMsIHRpbWVzID0gZHQsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmMgPSBTRUlSLCBwYXJtcyA9IHBhcm1zKSkKICAgICAgICAKICAgICAgICBob2xkZXIgPC0gcmJpbmQoaG9sZGVyLCBjYmluZCh0ID0gc2ltdWxhdGlvblssICJ0aW1lIl0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGluZmVjdGVkID0gc2ltdWxhdGlvblssICJJIl0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGFscGhhID0gYWxwaGEsIGdfcHJvcCA9IGdfcHJvcCkpCiAgICB9CiAgICAKfQpgYGAKCldlIG5vdyBoYXZlIGFuIG9iamVjdCB3aGljaCBjb250YWlucyB0aGUgaW5mZWN0aW91cyBjdXJ2ZSBmb3IgYWxsIDI1IGNvbWJpbmF0aW9ucyBvZiBpbmZlY3Rpb3VzIHBlcmlvZCBhbmQgaW5jdWJhdGlvbiBwZXJpb2RzLgoKV2UgY2FuIHBsb3QgaXQgbGlrZSBzbzoKCmBgYHtyfQpob2xkZXIgPC0gaG9sZGVyICU+JSAKICAgIGFzX3RpYmJsZSgpICU+JSAKICAgIG11dGF0ZShpbmZlY3Rpb3VzX2R1cmF0aW9uID0gcm91bmQoMS9hbHBoYSwgMCkpCgpnZ3Bsb3QoaG9sZGVyLCBhZXMoeCA9IHQsIHkgPSBpbmZlY3RlZCkpICsgCiAgICBnZW9tX2xpbmUoKSArIAogICAgZmFjZXRfZ3JpZChpbmZlY3Rpb3VzX2R1cmF0aW9uIH4gZ19wcm9wKSArIAogICAgc2NhbGVfeF9jb250aW51b3VzKGxpbWl0cyA9IGMoMCwgMTAwKSkgKyAKICAgIHRoZW1lX21pbmltYWwoKSArIAogICAgbGFicyh4ID0gIlRpbWUiLCB5ID0gIlByb3BvcnRpb24gSW5mZWN0ZWQiLCAKICAgICAgICAgdGl0bGUgPSAiSW5mZWN0aW91cyBjdXJ2ZXMgZm9yIHZhcmlvdXMgaW5mZWN0aW91cyBhbmQgaW5jdWJhdGlvbiBwZXJpb2RzIikKYGBgCgpFYWNoIHBhbmVsIHJlcHJlc2VudHMgc29tZSBpbnN0YW5jZSBvZiAkXGFscGhhJCBhbmQgJFxnYW1tYSQgZnJvbSBvdXIgc2V0cy4gRWFjaCByb3cgaXMgb25lIHNldCBvZiBpbmZlY3Rpb3VzIHBlcmlvZHMgKGkuZS4sIHRoZSB0b3AgaXMgYW4gaW5mZWNpdG91cyBwZXJpb2Qgb2YgMSBkYXksIHRoZSBib3R0b20gaXMgMjAgZGF5cykuIFRoZSBjb2x1bW5zIHNob3cgdGhlIGNoYW5naW5nIGluY3ViYXRpb24gcGVyaW9kcyBhcyBhIHJhdGlvcyBvZiB0aGUgaW5mZWN0aW91cyBwZXJpb2RzLiAKCklmIHlvdSBmb2N1cyBvbiBhbnkgcm93IGFuZCBnbyBmcm9tIGxlZnQgdG8gcmlnaHQsIHlvdSBjYW4gc2VlIHRoYXQgYXMgdGhlIGxhdGVudCBwZXJpb2QgZ2V0cyBsb25nZXIgKHJlbGF0aXZlIHRvIHRoZSBpbmZlY3Rpb3VzIHBlcmlvZCksIHRoZSBlcGlkZW1pYyBwZWFrIChtYXhpbXVtIG51bWJlciBvZiBpbmZlY3RlZCkgZ2V0cyBsb3dlciAqYW5kKiBhIGxpdHRsZSBtb3JlIGRlbGF5ZWQuIAoKSWYgeW91IGZvY3VzIG9uIGFueSBjb2x1bW4gYW5kIGdvIHRvcCB0byBib3R0b20sIHlvdSBjYW4gc2VlIHRoYXQgYXMgeW91IGluY3JlYXNlIHRoZSBkdXJhdGlvbiBvZiBpbmZlY3Rpb3VzbmVzcywgdGhlIGVwaWRlbWljIHBlYWsgaXMgc2hpZnRlZCB0byB0aGUgcmlnaHQgKGJlY2F1c2Ugd2UgYXJlIGFsc28gaW5jcmVhc2luZyB0aGUgYW1vdW50IG9mIHRpbWUgdGhleSBzdGF5IGluIHRoZSBsYXRlbnQgcGVyaW9kKSwgYW5kIHRoZSAicGVhayIgaXMgZ2VuZXJhbGx5IGZhdHRlciBiZWNhdXNlIHRoZSB0aW1lIHNwZW50IGluIHRoZSBjb21wYXJ0bWVudCBpcyBsb25nZXIuIEl0IGRvZXMgbm90LCBob3dldmVyLCBzZWVtIGxpa2UgdGhlcmUgaXMgYW55IGltcGFjdCBvbiB0aGUgaGVpZ2h0IG9mIHRoZSBwZWFrLiAKCiMjIFRha2UgaG9tZQoKTGF0ZW50IHBlcmlvZHMgbWF0dGVyIGFzIHRoZXkgYmVjb21lIGxvbmdlciByZWxhdGl2ZSB0byB5b3VyIGluZmVjdGlvdXMgcGVyaW9kLiBGb3IgbGF0ZW50IHBlcmlvZHMgbGFzdGluZyAkXGx0IDEwXCUkIG9mIHRoZSBpbmZlY3Rpb3VzIHBlcmlvZCwgbW9kZWxlcnMgd2lsbCBvZnRlbiBsZWF2ZSBpdCBvdXQgZm9yIHNpbXBsaWNpdHkuIFlvdSBzaG91bGQgYWxtb3N0IGNlcnRhaW5seSBpbmNsdWRlIGl0IGlmIGl0IGlzIGdyZWF0ZXIgdGhhbiB0aGF0IHRob3VnaC4KCgoK