Introduction

This is annotated code for the headache_ode.R file. For convenience, the code from that file is shown below (with comments removed):

library(deSolve)

times <- seq(0, 5000, by = 1)
yinit <- c(no_headache = 0.95, headache = 0.05)
parameters <- c(incidence = 0.02, recovery = 0.01, 
                birth = 0.1, death = 0.2)

headache_model <- function(times, yinit, parameters) {
    with(as.list(c(yinit, parameters)), {

        dno_headache <- recovery*headache - 
            incidence*no_headache - death*no_headache 

        dheadache <- incidence*no_headache - recovery*headache - 
            death*headache + birth*(headache+no_headache)

        comparts <- list(c(dno_headache, dheadache))

        return(comparts)
    })
}

result <- as.data.frame(ode(y = yinit, times = times, 
                            func = headache_model, parms = parameters))

matplot(x = result[, "time"], 
        y = result[, c("no_headache", "headache")], 
        type = "l")

Calling deSolve

library(deSolve)

To solve differential equations, we use the ode() command. This command comes from the deSolve package, which means we must load the package before we can run the command.

All scripts that use ode() will have the library(deSolve) command at the top.

Setting parameters

times <- seq(0, 5000, by = 1)
yinit <- c(no_headache = 0.95, headache = 0.05)
parameters <- c(incidence = 0.02, recovery = 0.01, 
                birth = 0.1, death = 0.2)

Like the discrete version, we need to set up initial values, parameter values, and a time sequence in order for ode() to work. Here we will use seq(0, 5000, by = 1) to specify that we want a vector of numbers from 0 to 5,000 in steps of 1.

print(yinit)
no_headache    headache 
       0.95        0.05 
print(parameters)
incidence  recovery     birth     death 
     0.02      0.01      0.10      0.20 

Creating the model

Suppose this is our model of headaches:

\[ \begin{align} \frac{d N}{dt} &= \gamma P - \beta N - dN \\ \frac{d P}{dt} &= \beta N - \gamma P - dP + b(N+P) \end{align} \] What is this model saying about births and deaths?

Unlike the discrete model where we created our model inside of a for() loop, for ODE models, we will create a function that contains our models. This function will require three inputs: (1) a time sequence, (2) initial values of each compartment, and (3) the values of any parameters. In addition, we’re going to call this function headache_model.

headache_model <- function(times, yinit, parameters) {
    with(as.list(c(yinit, parameters)), {
        dno_headache <- recovery*headache - 
            incidence*no_headache - death*no_headache 
        dheadache <- incidence*no_headache - recovery*headache - 
            death*headache + birth*(headache+no_headache)
        comparts <- list(c(dno_headache, dheadache))
        return(comparts)
    })
}

You will only be expected to modify things inside of the with block. Specifically the dno_headache, dheadache, comparts, and return() lines. The with() command simply allows us to refer to the named elements in our yinit and parameter variables without the need to specify the object they are in. For example, instead of specifying parameters$recovery * yinit$headache, we can simply say recovery*headache.

In order for ode() to work, we must return a list of compartments. We create comparts which stores our compartments as a list using the list() command.

Running our model

Now that we have our model (as a function), we can use ode() to solve it:

result <- as.data.frame(ode(y = yinit, times = times, 
                            func = headache_model, parms = parameters))

ode() will return a named list, but in general, we prefer data.frames which are easier to work with. Thus we simply run ode() and convert the output to a dataframe using as.data.frame(), storing the result in result.

Plot the results

Let’s see what ode() returned:

head(result)

Note that ode() appended a column named time to our model in addition to the compartments.

We’ll use matplot() again to plot the results:

matplot(x = result[, "time"], 
        y = result[, c("no_headache", "headache")], 
        type = "l")

Above, we referred to the columns by name, but we can also refer to them by position:

matplot(x = result[, 1], 
        y = result[, 2:3], 
        type = "l")

Trimming the x-axis

We can take a close look at the beginning of our model by specifying our x-axis limits with the xlim option. This option takes a vector of length 2 with the first element as the lower bound and the second element as the upper bound.

matplot(x = result[, 1], 
        y = result[, 2:3], 
        type = "l", 
        xlim = c(0, 50))

Try playing with other parameters.

LS0tCnRpdGxlOiAiQW5ub3RhdGVkIGNvZGUgZm9yIGBoZWFkYWNoZV9vZGUuUmAiCmF1dGhvcjogIk1hdGhldyBLaWFuZyIKZGF0ZTogImByIGZvcm1hdChTeXMudGltZSgpLCAnJUIgJWQsICVZJylgIgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgICB0b2M6IHRydWUKLS0tCgojIyBJbnRyb2R1Y3Rpb24KClRoaXMgaXMgYW5ub3RhdGVkIGNvZGUgZm9yIHRoZSBbYGhlYWRhY2hlX29kZS5SYF0oaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL21raWFuZy9lcGk1MDFfcmV2aWV3L21hc3Rlci9leGVyY2lzZXMvMDFfaGVhZGFjaGVzL2hlYWRhY2hlX29kZS5SKSBmaWxlLiBGb3IgY29udmVuaWVuY2UsIHRoZSBjb2RlIGZyb20gdGhhdCBmaWxlIGlzIHNob3duIGJlbG93ICh3aXRoIGNvbW1lbnRzIHJlbW92ZWQpOgoKYGBge3IgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShkZVNvbHZlKQoKdGltZXMgPC0gc2VxKDAsIDUwMDAsIGJ5ID0gMSkKeWluaXQgPC0gYyhub19oZWFkYWNoZSA9IDAuOTUsIGhlYWRhY2hlID0gMC4wNSkKcGFyYW1ldGVycyA8LSBjKGluY2lkZW5jZSA9IDAuMDIsIHJlY292ZXJ5ID0gMC4wMSwgCiAgICAgICAgICAgICAgICBiaXJ0aCA9IDAuMSwgZGVhdGggPSAwLjIpCgpoZWFkYWNoZV9tb2RlbCA8LSBmdW5jdGlvbih0aW1lcywgeWluaXQsIHBhcmFtZXRlcnMpIHsKICAgIHdpdGgoYXMubGlzdChjKHlpbml0LCBwYXJhbWV0ZXJzKSksIHsKCiAgICAgICAgZG5vX2hlYWRhY2hlIDwtIHJlY292ZXJ5KmhlYWRhY2hlIC0gCiAgICAgICAgICAgIGluY2lkZW5jZSpub19oZWFkYWNoZSAtIGRlYXRoKm5vX2hlYWRhY2hlIAoKICAgICAgICBkaGVhZGFjaGUgPC0gaW5jaWRlbmNlKm5vX2hlYWRhY2hlIC0gcmVjb3ZlcnkqaGVhZGFjaGUgLSAKICAgICAgICAgICAgZGVhdGgqaGVhZGFjaGUgKyBiaXJ0aCooaGVhZGFjaGUrbm9faGVhZGFjaGUpCgogICAgICAgIGNvbXBhcnRzIDwtIGxpc3QoYyhkbm9faGVhZGFjaGUsIGRoZWFkYWNoZSkpCgogICAgICAgIHJldHVybihjb21wYXJ0cykKICAgIH0pCn0KCnJlc3VsdCA8LSBhcy5kYXRhLmZyYW1lKG9kZSh5ID0geWluaXQsIHRpbWVzID0gdGltZXMsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuYyA9IGhlYWRhY2hlX21vZGVsLCBwYXJtcyA9IHBhcmFtZXRlcnMpKQoKbWF0cGxvdCh4ID0gcmVzdWx0WywgInRpbWUiXSwgCiAgICAgICAgeSA9IHJlc3VsdFssIGMoIm5vX2hlYWRhY2hlIiwgImhlYWRhY2hlIildLCAKICAgICAgICB0eXBlID0gImwiKQpgYGAKCiMjIENhbGxpbmcgYGRlU29sdmVgCmBgYHtyfQpsaWJyYXJ5KGRlU29sdmUpCmBgYAoKVG8gc29sdmUgZGlmZmVyZW50aWFsIGVxdWF0aW9ucywgd2UgdXNlIHRoZSBgb2RlKClgIGNvbW1hbmQuIFRoaXMgY29tbWFuZCBjb21lcyBmcm9tIHRoZSBgZGVTb2x2ZWAgcGFja2FnZSwgd2hpY2ggbWVhbnMgd2UgbXVzdCBsb2FkIHRoZSBwYWNrYWdlIGJlZm9yZSB3ZSBjYW4gcnVuIHRoZSBjb21tYW5kLiAKCkFsbCBzY3JpcHRzIHRoYXQgdXNlIGBvZGUoKWAgd2lsbCBoYXZlIHRoZSBgbGlicmFyeShkZVNvbHZlKWAgY29tbWFuZCBhdCB0aGUgdG9wLgoKIyMgU2V0dGluZyBwYXJhbWV0ZXJzCmBgYHtyfQp0aW1lcyA8LSBzZXEoMCwgNTAwMCwgYnkgPSAxKQp5aW5pdCA8LSBjKG5vX2hlYWRhY2hlID0gMC45NSwgaGVhZGFjaGUgPSAwLjA1KQpwYXJhbWV0ZXJzIDwtIGMoaW5jaWRlbmNlID0gMC4wMiwgcmVjb3ZlcnkgPSAwLjAxLCAKICAgICAgICAgICAgICAgIGJpcnRoID0gMC4xLCBkZWF0aCA9IDAuMikKYGBgCgpMaWtlIHRoZSBkaXNjcmV0ZSB2ZXJzaW9uLCB3ZSBuZWVkIHRvIHNldCB1cCBpbml0aWFsIHZhbHVlcywgcGFyYW1ldGVyIHZhbHVlcywgYW5kIGEgdGltZSBzZXF1ZW5jZSBpbiBvcmRlciBmb3IgYG9kZSgpYCB0byB3b3JrLiBIZXJlIHdlIHdpbGwgdXNlIGBzZXEoMCwgNTAwMCwgYnkgPSAxKWAgdG8gc3BlY2lmeSB0aGF0IHdlIHdhbnQgYSB2ZWN0b3Igb2YgbnVtYmVycyBmcm9tIDAgdG8gNSwwMDAgaW4gc3RlcHMgb2YgMS4gCgotIFVubGlrZSB0aGUgZGlzY3JldGUgbW9kZWwsIHdlIHdvbid0IG5lZWQgdG8gdXNlIGludGVnZXJzIGFuZCB3ZSBjb3VsZCBoYXZlIGRvbmUgYHNlcSgwLCA1MDAwLCBieSA9IC4wMDAxKWAuIFRoaXMgd291bGQgcmVzdWx0IGluIG91ciBlcGlkZW1pYyBjdXJ2ZXMgYmVpbmcgbXVjaCBzbW9vdGhlciBiZWNhdXNlIGBvZGUoKWAgd291bGQgZXZhbHVhdGUgb3VyIG1vZGVsIGF0IGV2ZXJ5IC4wMDAxIHVuaXRzIGluc3RlYWQgb2YgZXZlcnkgMSB1bml0LiBIb3dldmVyLCB0aGlzIHdvdWxkIGFsc28gcmVzdWx0IGluIG11Y2ggbG9uZ2VyIGNvbXB1dGF0aW9uIHRpbWUuIENob29zaW5nIHRoZSBsZW5ndGggb2YgeW91ciB0aW1lIHNlcXVlbmNlIGFuZCB0aGUgcmVzb2x1dGlvbiB3aWxsIGJlIGEgYmFsYW5jaW5nIGFjdCB5b3UnbGwgbmVlZCB0byBmaWd1cmUgb3V0IHRocm91Z2hvdXQgdGhlIGNvdXJzZSwgYnV0IHdlIHdpbGwgZ2VuZXJhbGx5IGdpdmUgeW91IGdvb2Qgc3RhcnRpbmcgdmFsdWVzLgoKLSBXZSBzdG9yZSBhbGwgb2Ygb3VyIGluaXRpYWwgdmFsdWVzIGluIGEgdmFyaWFibGUgY2FsbGVkIGB5aW5pdGAuIFRoaXMgaXMgYSBuYW1lZCB2ZWN0b3Igc28gdGhhdCB3aGVuIHdlIGNhbGwgYHlpbml0YCBpdCB3aWxsIHJlc3VsdCBpbiBvbmUgZWxlbWVudCBjYWxsZWQgYG5vX2hlYWRhY2hlYCBhbmQgb25lIGNhbGxlZCBgaGVhZGFjaGVgLiBXaGVuIHlvdSBhZGQgb3IgcmVtb3ZlIGNvbXBhcnRtZW50cywgeW91J2xsIG5lZWQgdG8gYWRqdXN0IHRoZSBgeWluaXRgIHZhcmlhYmxlIGFjY29yZGluZ2x5LiBUaGVyZSBzaG91bGQgYWx3YXlzIGJlIGFzIG1hbnkgZWxlbWVudHMgYXMgdGhlcmUgYXJlIGNvbXBhcnRtZW50cy4KCmBgYHtyfQpwcmludCh5aW5pdCkKYGBgCgotIFNpbWlsYXJseSwgd2Ugd2lsbCBjcmVhdGUgYSBuYW1lZCB2ZWN0b3IgdG8gc3RvcmUgYWxsIG91ciBwYXJhbWV0ZXJzLiBXaGVuIHlvdSBtb2RpZnkgeW91ciBjb21wYXJ0bWVudHMgd2l0aCBuZXcgcGFyYW1ldGVycywgeW91J2xsIG5lZWQgdG8gc3BlY2lmeSB0aGUgdmFsdWUgb2YgdGhlc2UgcGFyYW1ldGVycyBpbiB0aGlzIHZhcmlhYmxlLgoKYGBge3J9CnByaW50KHBhcmFtZXRlcnMpCmBgYAoKIyMgQ3JlYXRpbmcgdGhlIG1vZGVsCgpTdXBwb3NlIHRoaXMgaXMgb3VyIG1vZGVsIG9mIGhlYWRhY2hlczoKCiQkClxiZWdpbnthbGlnbn0KXGZyYWN7ZCBOfXtkdH0gJj0gXGdhbW1hIFAgLSBcYmV0YSBOIC0gZE4gXFwKXGZyYWN7ZCBQfXtkdH0gJj0gXGJldGEgTiAtIFxnYW1tYSBQIC0gZFAgKyBiKE4rUCkKXGVuZHthbGlnbn0KJCQKV2hhdCBpcyB0aGlzIG1vZGVsIHNheWluZyBhYm91dCBiaXJ0aHMgYW5kIGRlYXRocz8KClVubGlrZSB0aGUgZGlzY3JldGUgbW9kZWwgd2hlcmUgd2UgY3JlYXRlZCBvdXIgbW9kZWwgaW5zaWRlIG9mIGEgYGZvcigpYCBsb29wLCBmb3IgT0RFIG1vZGVscywgd2Ugd2lsbCBjcmVhdGUgYSBmdW5jdGlvbiB0aGF0IGNvbnRhaW5zIG91ciBtb2RlbHMuIFRoaXMgZnVuY3Rpb24gd2lsbCByZXF1aXJlIHRocmVlIGlucHV0czogKDEpIGEgdGltZSBzZXF1ZW5jZSwgKDIpIGluaXRpYWwgdmFsdWVzIG9mIGVhY2ggY29tcGFydG1lbnQsIGFuZCAoMykgdGhlIHZhbHVlcyBvZiBhbnkgcGFyYW1ldGVycy4gSW4gYWRkaXRpb24sIHdlJ3JlIGdvaW5nIHRvIGNhbGwgdGhpcyBmdW5jdGlvbiBgaGVhZGFjaGVfbW9kZWxgLgoKYGBge3J9CmhlYWRhY2hlX21vZGVsIDwtIGZ1bmN0aW9uKHRpbWVzLCB5aW5pdCwgcGFyYW1ldGVycykgewogICAgd2l0aChhcy5saXN0KGMoeWluaXQsIHBhcmFtZXRlcnMpKSwgewoKICAgICAgICBkbm9faGVhZGFjaGUgPC0gcmVjb3ZlcnkqaGVhZGFjaGUgLSAKICAgICAgICAgICAgaW5jaWRlbmNlKm5vX2hlYWRhY2hlIC0gZGVhdGgqbm9faGVhZGFjaGUgCgogICAgICAgIGRoZWFkYWNoZSA8LSBpbmNpZGVuY2Uqbm9faGVhZGFjaGUgLSByZWNvdmVyeSpoZWFkYWNoZSAtIAogICAgICAgICAgICBkZWF0aCpoZWFkYWNoZSArIGJpcnRoKihoZWFkYWNoZStub19oZWFkYWNoZSkKCiAgICAgICAgY29tcGFydHMgPC0gbGlzdChjKGRub19oZWFkYWNoZSwgZGhlYWRhY2hlKSkKCiAgICAgICAgcmV0dXJuKGNvbXBhcnRzKQogICAgfSkKfQpgYGAKCllvdSB3aWxsIG9ubHkgYmUgZXhwZWN0ZWQgdG8gbW9kaWZ5IHRoaW5ncyBpbnNpZGUgb2YgdGhlIGB3aXRoYCBibG9jay4gU3BlY2lmaWNhbGx5IHRoZSBgZG5vX2hlYWRhY2hlYCwgYGRoZWFkYWNoZWAsIGBjb21wYXJ0c2AsIGFuZCBgcmV0dXJuKClgIGxpbmVzLiBUaGUgYHdpdGgoKWAgY29tbWFuZCBzaW1wbHkgYWxsb3dzIHVzIHRvIHJlZmVyIHRvIHRoZSBuYW1lZCBlbGVtZW50cyBpbiBvdXIgYHlpbml0YCBhbmQgYHBhcmFtZXRlcmAgdmFyaWFibGVzIHdpdGhvdXQgdGhlIG5lZWQgdG8gc3BlY2lmeSB0aGUgb2JqZWN0IHRoZXkgYXJlIGluLiBGb3IgZXhhbXBsZSwgaW5zdGVhZCBvZiBzcGVjaWZ5aW5nIGBwYXJhbWV0ZXJzJHJlY292ZXJ5ICogeWluaXQkaGVhZGFjaGVgLCB3ZSBjYW4gc2ltcGx5IHNheSBgcmVjb3ZlcnkqaGVhZGFjaGVgLgoKSW4gb3JkZXIgZm9yIGBvZGUoKWAgdG8gd29yaywgd2UgbXVzdCBgcmV0dXJuYCBhIGxpc3Qgb2YgY29tcGFydG1lbnRzLiBXZSBjcmVhdGUgYGNvbXBhcnRzYCB3aGljaCBzdG9yZXMgb3VyIGNvbXBhcnRtZW50cyBhcyBhIGxpc3QgdXNpbmcgdGhlIGBsaXN0KClgIGNvbW1hbmQuCgojIyBSdW5uaW5nIG91ciBtb2RlbApOb3cgdGhhdCB3ZSBoYXZlIG91ciBtb2RlbCAoYXMgYSBmdW5jdGlvbiksIHdlIGNhbiB1c2UgYG9kZSgpYCB0byBzb2x2ZSBpdDoKCmBgYHtyfQpyZXN1bHQgPC0gYXMuZGF0YS5mcmFtZShvZGUoeSA9IHlpbml0LCB0aW1lcyA9IHRpbWVzLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmMgPSBoZWFkYWNoZV9tb2RlbCwgcGFybXMgPSBwYXJhbWV0ZXJzKSkKCmBgYAoKYG9kZSgpYCB3aWxsIHJldHVybiBhIG5hbWVkIGxpc3QsIGJ1dCBpbiBnZW5lcmFsLCB3ZSBwcmVmZXIgYGRhdGEuZnJhbWVgcyB3aGljaCBhcmUgZWFzaWVyIHRvIHdvcmsgd2l0aC4gVGh1cyB3ZSBzaW1wbHkgcnVuIGBvZGUoKWAgYW5kIGNvbnZlcnQgdGhlIG91dHB1dCB0byBhIGRhdGFmcmFtZSB1c2luZyBgYXMuZGF0YS5mcmFtZSgpYCwgc3RvcmluZyB0aGUgcmVzdWx0IGluIGByZXN1bHRgLgoKCiMjIFBsb3QgdGhlIHJlc3VsdHMKTGV0J3Mgc2VlIHdoYXQgYG9kZSgpYCByZXR1cm5lZDoKYGBge3J9CmhlYWQocmVzdWx0KQpgYGAKCk5vdGUgdGhhdCBgb2RlKClgIGFwcGVuZGVkIGEgY29sdW1uIG5hbWVkIGB0aW1lYCB0byBvdXIgbW9kZWwgaW4gYWRkaXRpb24gdG8gdGhlIGNvbXBhcnRtZW50cy4KCldlJ2xsIHVzZSBgbWF0cGxvdCgpYCBhZ2FpbiB0byBwbG90IHRoZSByZXN1bHRzOgoKYGBge3J9Cm1hdHBsb3QoeCA9IHJlc3VsdFssICJ0aW1lIl0sIAogICAgICAgIHkgPSByZXN1bHRbLCBjKCJub19oZWFkYWNoZSIsICJoZWFkYWNoZSIpXSwgCiAgICAgICAgdHlwZSA9ICJsIikKYGBgCgpBYm92ZSwgd2UgcmVmZXJyZWQgdG8gdGhlIGNvbHVtbnMgYnkgbmFtZSwgYnV0IHdlIGNhbiBhbHNvIHJlZmVyIHRvIHRoZW0gYnkgcG9zaXRpb246CmBgYHtyLCBldmFsPUZBTFNFfQptYXRwbG90KHggPSByZXN1bHRbLCAxXSwgCiAgICAgICAgeSA9IHJlc3VsdFssIDI6M10sIAogICAgICAgIHR5cGUgPSAibCIpCmBgYAoKIyMjIFRyaW1taW5nIHRoZSB4LWF4aXMKV2UgY2FuIHRha2UgYSBjbG9zZSBsb29rIGF0IHRoZSBiZWdpbm5pbmcgb2Ygb3VyIG1vZGVsIGJ5IHNwZWNpZnlpbmcgb3VyIHgtYXhpcyBsaW1pdHMgd2l0aCB0aGUgYHhsaW1gIG9wdGlvbi4gVGhpcyBvcHRpb24gdGFrZXMgYSB2ZWN0b3Igb2YgbGVuZ3RoIDIgd2l0aCB0aGUgZmlyc3QgZWxlbWVudCBhcyB0aGUgbG93ZXIgYm91bmQgYW5kIHRoZSBzZWNvbmQgZWxlbWVudCBhcyB0aGUgdXBwZXIgYm91bmQuCgpgYGB7cn0KbWF0cGxvdCh4ID0gcmVzdWx0WywgMV0sIAogICAgICAgIHkgPSByZXN1bHRbLCAyOjNdLCAKICAgICAgICB0eXBlID0gImwiLCAKICAgICAgICB4bGltID0gYygwLCA1MCkpCmBgYAoKVHJ5IHBsYXlpbmcgd2l0aCBvdGhlciBwYXJhbWV0ZXJzLg==