Open presentation in your browser

Why?

• Access to links, downloads, source code

How?

• Find my home page via web search or navigate to https://wwwpub.zih.tu-dresden.de/~kneis
• Find “Hacking limnology” in the “teaching” section
• Click “open” and follow the 2^nd link

Basic software

• Statistical computing software “R”

• RStudio or a good text editor like geany

• Spreadsheet software to edit workbooks in “.xlsx” format (e.g., LibreOffice “calc”)

• gfortran compiler (part of the RTools for Windows)

Is gfortran installed?

Add-on packages for R

• We will primarily use one package: rodeo

• Dependencies on other packages (e.g., deSolve) are resolved automatically upon installation

• We will install the latest version of rodeo in a moment

Copying code

• Select / create a working directory

• Therein, create a text file ending with “.r” or “.R”)

• Copy & paste R code in given order

# to copy this code to the clipboard, place
# the mouse over the box & click the icon in
# the top-right corner
print("I'm a piece of R code. Copy me!")

Executing code

RStudio

• Have file with R code open

• Session \(\rightarrow\) Set work. dir. \(\rightarrow\) To source file loc.

• Hit source button or press Ctrl + Shift + S

Using stop() and print()

• Do NOT run scripts line by line but execute all statements from top of the file!

• Inspect intermediate results by print() & stop()

myData <- runif(999)

print(head(myData))
stop("let me check if data is good")

myOutputs <- sample(myData, size=5)

To disable the stop(), just put a # in front.

Objectives

• Simplest approach to a rodeo-based model

• Learn how to extend an existing model

• Learn how to update model inputs

• Understand the use of dynamic forcings

Continuous transfer culture

Equations

\[ \begin{align} \tfrac{d}{dt} B & = & & \color{teal}{\mu \cdot B \cdot \left( \tfrac{R}{R+h} \right)} & \color{blue}{- \tfrac{1}{\tau} \cdot B} \\ \tfrac{d}{dt} R & = & \color{red}{-\tfrac{1}{y}} \cdot & \color{teal}{\mu \cdot B \cdot \left( \tfrac{R}{R+h} \right)} & \color{blue}{+ \tfrac{1}{\tau} \cdot (R_{in} - R)} \end{align} \]

Representation as tables

• Choose / create a working directory

• Download the workbook reactor_eqns_v1.xlsx into that directory

• Inspect contents using spreadsheet software

Declarations (variables)

Declarations (parameters)

Equations

Overview of steps

• Build model object

• Set / check input data

• Perform integration

• Inspect results

Building the model

1. Go to the folder with the downloaded workbook

2. Create & open a new script file there (e.g. “reactor.r”)

3. Adjust R’s working directory (if needed)

rm(list=ls())
library("rodeo")

m <- buildFromWorkbook("reactor_eqns_v1.xlsx")
print(class(m))

[1] "rodeo" "R6"   

OO-syntax

Usual R code

y <- fun( x, args )    # don't run !

• x is a variable passed to function fun

Using R6 classes

y <- x$fun( args )     # don't run !

• x is an object having a member function fun

Check attached data

print("Current initial values:")
print( m$getVars() )

[1] "Current initial values:"
     B      R 
100000     21 

print("Current values of parameters:")
print( m$getPars() )

[1] "Current values of parameters:"
    tau      mu       y       h     Rin 
2.4e+01 1.0e+00 1.0e+07 2.1e+00 2.1e+01 

Running a simulation

x <- m$dynamics(times=0:24, fortran=F)

print(head(x,3))
print(tail(x,3))

     time        B        R    growth    renewal
[1,]    0 100000.0 21.00000  90909.09 0.04166667
[2,]    1 238071.7 20.98578 216415.47 0.04166667
[3,]    2 566731.0 20.95253 515103.91 0.04166667
      time         B          R  growth    renewal
[23,]   22 209127187 0.09128602 8711957 0.04166667
[24,]   23 209125546 0.09128677 8711957 0.04166667
[25,]   24 209123971 0.09128748 8711957 0.04166667

Visualization

plt <- function(x, m) { # x: data, m: model obj.
  
  vn <- m$namesVars()   # variable names
  nv <- m$lenVars()     # number of variables
  
  matplot(x[,"time"], x[,vn], log="y", type="l",
    lty=1:nv, col=1:nv, xlab="Hour", ylab="Value")
  legend("right", bty="n", lty=1:nv,
    col=1:nv, legend=vn)
}

plt(x, m)

Competition

Modifications needed

• Add drug resistant bacteria as new state variable \(BR\)

• Add antibiotic as new state variable \(A\)

• Account for input of antibiotic (\(Ain\)) & removal of \(A\)

• Account for growth of \(BR\) and the “cost” of resistance implying \(\mu(BR) < \mu(B)\)

• Account for drug susceptibility of strain \(B\)

• Consider removal of \(BR\) during transfer

Modifying the workbook

Usual procedure

• Make a copy of “reactor_eqns_v1.xlsx”

• Save as “reactor_eqns_v2.xlsx”

• Implement the modifications

Fallback

• Download the ready-made version of the workbook reactor_eqns_v2.xlsx

Declarations (variables)

Declarations (parameters)

Equations

Re-build & simulate

m <- buildFromWorkbook("reactor_eqns_v2.xlsx")

x <- m$dynamics(times=0:48, fortran=F)

print(x[1:3, 1:6])

     time        B       BR        R A  growth_B
[1,]    0 100000.0 100000.0 21.00000 0  90909.09
[2,]    1 238066.7 217379.7 20.97364 0 216399.53
[3,]    2 566666.6 472471.2 20.91449 0 514960.05

Visualization

We reuse the previously defined function

plt(x, m)

Changing model inputs

• Option 1: Modify defaults in worksheet

• Option 2: Modify through rodeo methods

p <- m$getPars()    # current settings
p["Ain"] <- 10      # now with antibiotic
m$setPars(p)        # update settings

# no need to rebuild/recompile the model
x <- m$dynamics(times=0:48, fortran=F)
plt(x, m)

Example: Varying drug exposure

1. Declare Ain as a function, not as a parameter

2. Replace Ain by Ain(time) in expressions (time is a reserved word)

3. Implement the function in R (or Fortran)

Modifying the workbook

Usual procedure

• Make a copy of “reactor_eqns_v2.xlsx”

• Save as “reactor_eqns_v3.xlsx”

• Modify as necessary

Fallback

• Download the ready-made version of the workbook reactor_eqns_v3.xlsx

Declarations (functions)

Declarations (parameters)

Equations

Function implementation in R

Ain <- function(time) {
  h <- floor((time/24 - floor(time/24)) * 24)
  if (h %in% c(0:5))
    10
  else
    0
}

\(\rightarrow\) Antibiotic added from midnight till 6 am

Re-build & simulate

m <- buildFromWorkbook("reactor_eqns_v3.xlsx")

x <- m$dynamics(times=0:240, fortran=F)

plt(x, m)

Objectives

• Solution of PDE problems with ODE strategies

• Use of Fortran for improved performance

Streeter & Phelps-like model

This particular version

• Allows for multiple sources & time-varying inputs

• Treats bacteria as a dynamic state variable

Reactive transport is a PDE problem

For lateral & vertical mixing, a 1D model may fit

\[ \frac{\partial c}{\partial t} = \color{teal}{-u \cdot \frac{\partial c}{\partial x}} + \color{blue}{\frac{\partial}{\partial x} \left( D \cdot \frac{\partial c}{\partial x} \right)} + \color{red}{r(c, p)} + \color{grey}{s} \]

Method of Lines (MOL)

Method of Lines (MOL)

Method of Lines (MOL)

\[\frac{\partial c_{\color{red}i}}{\partial t} = -u \cdot \frac{c_{\color{red}{i}} - c_{\color{blue}{i-1}}}{x_{\color{red}{i}} - x_{\color{blue}{i-1}}}\]

• Semi-discretization transforms PDE into a set of ODE

• The ODE for any cell \(\color{red}{i}\) references the state of neighboring cells (here: \(\color{blue}{i-1}\))

• ODE are simultaneous across cells (at least)

Method of Lines (MOL)

Finite difference schemes

• Backward, forward, central schemes

• Schemes differ w.r.t. stability, accuracy, and numerical dispersion

• See the documentation of function ‘fiadeiro’ in the ReacTran package

Representation as tables

• Choose / create a working directory

• Download the workbook river_eqns_v1.xlsx into that directory

• Inspect contents using spreadsheet software

Note that the ODE are specified for a single cell only, no matter how many cells the model consists of.

Declaration (variables)

Declarations (biological param.)

Declarations (hydraulic param.)

Declarations (hydraulic param.)

Declarations (boundary cond. param.)

Declarations (mask param.)

Equations

Equations

Overview of steps

• Build model object

• Set input data

• Perform integration

• Inspect results

Building the model

library("rodeo")

ncells <- 100     # number of cells
fortran <- TRUE   # TRUE requires compiler !

# import equations and declarations,
# translate into source code, and compile
m <- buildFromWorkbook("river_eqns_v1.xlsx",
  dim=ncells, fortran=fortran, na="NA")

Setting input data

Preparation of parameters

# retrieve defaults from workbook
p <- m$getParsTable()
# save as vector 
p <- setNames(p$default, p$name)
# copy to all cells (vector -> matrix)
p <- sapply(p, rep, length.out=ncells)

print(p[1:2, 1:5])

      mu   hs       hx   f   q
[1,] 0.1 0.01 0.015625 0.1 100
[2,] 0.1 0.01 0.015625 0.1 100

Setting input data

Preparation of parameters (cont.)

# adjust mask parameters
p[,"is_upstr"] <- c(1, rep(0, ncells-1))
p[,"is_centr"] <- c(0, rep(1, ncells-1))
p[,"has_src"] <- rep(0, ncells)
p[round(ncells) / 10, "has_src"] <- 1

# assign parameters to model object
m$setPars(p)

Setting input data

Similar procedure for initial values

v <- m$getVarsTable()
v <- setNames(v$default, v$name)
v <- sapply(v, rep, length.out=ncells)

m$setVars(v)

Running the simulation

times <- c(0, 12, 24, 48, 96, 120)
x <- m$dynamics(times=times, fortran=fortran)

print(x[ ,1:7])

     time   B.1   B.2   B.3   B.4   B.5   B.6
[1,]    0 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04
[2,]   12 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04
[3,]   24 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04
[4,]   48 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04
[5,]   96 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04
[6,]  120 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04

Transform results into 3D array

# dimension 1: time
# dimension 2: cell
# dimension 3: variables and process rates
a <- array(x[,colnames(x) != "time"],
  dim=c(nrow(x), ncells, sum(m$lenVars(),
    m$lenPros())),
  dimnames=list(time=x[,"time"], cell=1:ncells,
    variable=c(m$namesVars(), m$namesPros()))
)

print(a["48", "10", "X"])

[1] 0.3122861

Longitudinal profiles

plot_longit <- function(v, a, leg) {
  nc <- dim(a)[names(dimnames(a)) == "cell"]
  plot(c(1, nc), range(a[,,v]),
    type="n", xlab="Cell", ylab="")
  times <- dimnames(a)[["time"]]
  for (t in times)
    lines(1:nc, a[t, , v], col=match(t, times))
  if (leg) legend("right", bty="n", title="Time",
    lty=1, col=1:length(times), legend=times)
  abline(v=which(p[,"has_src"] == 1), lty=2)
  mtext(side=3, paste("Var.:",v), cex=par("cex"))
}

Longitudinal profiles

layout(matrix(1:3, nrow=1))

for (v in c("S", "X", "B")) {
  plot_longit(v, a, leg=(v=="S"))
}

layout(1)

Motivation

• ODEs may involve complicated expressions

• These can be wrapped into functions for clarity

• Let’s try with the Monod term …

Steps required

1. Declare the new function

2. Use the new function in equations

3. Implement the function (now in Fortran)

4. Rebuild / recompile the model

Step 1

• Make a copy of the workbook “river_eqns_v1.xlsx” and save as “river_eqns_v2.xlsx”

• On sheet “declarations”, declare the new function:

Step 2

On worksheet “equations” …

Step 3

Create a file “river_func.f95” and enter the following piece of Fortran

module functions
implicit none
contains

double precision function monod (x, h)
  double precision, intent(in):: x, h
  monod = x / (x + h)
end function

end module

Step 4

Rebuild / compile the model object using …

• the name of the modified workbook

• the “sources” argument to pass Fortran code

m <- buildFromWorkbook(
  "river_eqns_v2.xlsx",
  dim=ncells, fortran=fortran, na="NA",
  sources="river_func.f95"
)

In that order

• Use shown examples as guidance

• Look at the package vignette for more examples

• Provide me with a fully self-contained minimum working example, including …

  • Brief description of problem

  • Worksheets as tab-delimited text files (not .xlsx)

  • R / Fortran code